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PREFACE 
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Air  Force  Geophysics  Laboratory 
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Bedford,  Massachusetts  01731 

The  writers  express  their  thanks  to  Ms.  Eunice 
C.  Cronin,  Branch  Chief,  and  to  Mr.  John  F.  Kellaher  and  to  Mr. 
Austin  A.  Almon,  Contract  Monitors,  whose  support  and  tech- 
nical guidance  were  invaluable  in  the  performing  of  the  tasks 
described  in  this  report. 

In  addition,  the  authors  would  like  to  thank  the 
AFGL  research  personnel  for  whom  the  tasks  were  performed. 

Their  cooperation  and  technical  assistance  is  truly  appre- 
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1,  Introduction 


The  purpose  of  this  report  is  to  describe  most  of  the  mathematical 
analysis,  data  analysis  and  computer  programming  problems  performed  under 
Contract  Number  F19628-76-C-0203 . The  problems  vary  in  complexity  from 
straightforward  program  adaptation  to  tasks  requiring  analysis,  determining, 
implementing,  and  sometimes  deriving  algorithm  best  suited  to  perform  the 
calculations.  The  problems  discussed  are  done  so  in  summary  form.  The 
analysis  and  programming  techniques  required  are  outlined.  It  should  be 
noted  that  some  of  these  problems  are  still  being  analyzed  because  of  their 
complexity  or  continuing  nature.  Others  are  active  because  they  were  be- 
gun shortly  before  the  writing  of  this  report. 

In  the  subject  description,  problem  numbers  are  sometimes  suppli- 
mented  with  an  A or  B.  This  is  to  indicate  multiple  requests  under  the 
same  problem  number.  Problem  descriptions  which  are  for  follow-on  tasks 
or  slight  deviations  from  another  task  are  described  immediately  after  the 
primary  task.  Total  request  description  is  not  always  given,  rather  only 
the  modifications  are  discussed. 

Software  developed  for  tasks  completed  have  been  documented 
separately.  These  programs  can  be  obtained  from  The  Analysis  and  Simu- 
lation Branch  (SUA) , Air  Force  Geophysics  Laboratory  (AFGL)  upon  request 
by  referencing  the  appropriate  problem  numbers  and  project  numbers  listed 
with  each  task  description. 
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2.  Probabilities  and  Correlations  prom  Radar  Data 


INITIATOR:  I.  GRINGORTEN 


PROBLEM  NO.:  4862A  PROJECT  NO. : 8624 


BACKGROUND 

The  purpose  of  this  task  was  to  statistically  analyze  the  areal 
coverage,  determined  by  radar  echoes,  of  eight  equal  areas.  The  radar  echoes 
were  caused  by  precipitation.  The  radar  pictures  give  the  horizontal  dis- 
tribution of  precipitation  within  a given  radius  from  the  radar  station. 

There  is  a set  of  data  from  five  stations  for  each  day  of  the  four  mid- 
season months  (January,  April,  July  and  October)  in  1969-1974,  for  each 
3-hourly  interval.  The  data  gives  the  amount  of  areal  coverage  by  radar 
echoes  of  64  cells  of  equal  size  surrounding  a radar  station. 

With  the  64  cells  assembled  into  eight  groups,  the  task  was  to 
find  the  frequency  distribution  and  the  cumulative  probability  distribution 
of  areal  coverage  in  each  group.  Then  using  normalized  probabilities,  the 
correlation  between  the  groups  was  determined.  In  addition  to  the  standard 
correlation,  lag  correlations  between  the  seven  groups  used  as  predictors 
(groups  2 to  8)  and  the  first  group  were  calculated. 


ANALYSIS  AND  RESULTS 

Data  on  tape  was  available  from  five  different  stations,  namely 


STATION 

1 

2 

3 

4 

5 


NAME 

Minneapolis,  MN 
Key  West,  FL 
Wichita,  KS 
Cape  Hatteras,  NC 
Evansville,  IN 


2 


At  each  station,  radar  PPI  pictures  were  taken  once  every  three  hours  show- 
ing the  horizontal  distribution  of  rain  clouds  around  the  radar  station.  A 
copy  of  a transparent  template  (Figure  2-1)  was  placed  over  each  radar  pic- 
ture and  the  projection  was  adjusted  so  that  the  outermost  circle  fitted  on 
the  125-nmi  circle  of  the  template.  The  second  largest  circle  of  the  template 
represents  100-nmi,  the  innermost  circle  has  a scale-radius  of  25-nmi.  Be- 
tween the  25-nmi  and  the  100-nmi  circles,  the  area  is  divided  by  the  remaining 
three  circles  and  the  radial  lines  into  64  equal  cells  numbered  from  11  to 
74.  The  cells  numbered  11-18  form  Group  1;  those  numbered  19-26  form  Group 
2;  from  27-37,  Group  3;  and  so  on,  for  a total  of  eight  groups. 

In  each  radar  picture,  for  each  of  the  three  hours  throughout  the 
day,  for  each  day  of  the  four  months,  for  the  six  years,  the  fractional 
coverage  of  each  cell  (11-74)  was  determined  to  be  either  no  coverage,  or 
less  than  1/16,  1/8,  2/8,  up  to  8/8  (full  coverage) , and  recorded  respective- 
ly as  a blank,  0,  1,  2,  up  to  8 in  that  column  of  a punch  card  having  the 
same  number  as  the  cell  number. 

Software  was  written  to  combine  the  cells  into  the  eight  groups. 

Each  group  could  have  0 to  64  octas  of  rain-cloud  coverage.  (One  octa  is 
defined  as  1/8  coverage.)  The  frequency  distribution  for-  each  group  was  found 
in  terms  of  the  value  of  the  cells.  Then  the  cumulative  probability  that 
the  number  of  octas  exceeds  64,  63,  ...,  0 was  calculated.  The  case  where 
all  eight  cells  in  a group  were  blank  was  handled  separately.  After  the 
areal  coverage  in  Groups  1 to  8 was  determined  for  each  radar  picture,  the 
correlation  between  groups  was  determined.  Then  the  lag  correlations  between 
Group  1 and  the  eight  qroups  were  determined.  The  laqs  of  interest  were 
three-,  six-,  nine-  and  twelve-hour  lags. 

The  input  data  for  each  station  was  on  magnetic  tape.  The  informa- 
tion for  each  of  three  hours  for  every  day  in  the  four  months  for  the  six 
years  was  provided  on  the  tape. 


The  software  written  provided  the  following  information  (in  tabular 


output)  for  each  season/station  combination: 


Table  1 The  cumulative  relative  frequency  of  coverage  equal 

to  or  greater  than  L OCTAS  in  each  group,  for  L = 64, 
63 1,  0. 


Table  2 
Table  3 
Table  4 
Table  5 
Table  6 


the  cumulative  probability  distribution  for  Table  (1) . 
The  normalized  probability  values  of  Table  (2) . 

The  matrix  of  the  correlation  coefficients  (o. .). 

n 

The  matrix  of  the  lag  correlation  coefficients. 

Number  of  points  used  to  calculate  the  various 
correlations . 


DATA  PREPARATION 

In  preparing  the  data  to  be  run  through  the  program,  errors  in  the 
data  base  were  detected.  These  errors  were  corrected  before  the  program 
was  run.  The  following  is  a description  of  the  procedure  used  to  prepare 
the  total  data  base. 

The  data  for  five  stations  was  originally  on  9-track  and  7-track 
tapes.  The  first  step  was  to  convert  the  9-track  tapes  into  7-track  tapes. 
All  stations  contained  missing  and/or  mispunched  data.  The  corrections  were 
made  by  either  merging  two  records,  by  deleting  records,  or  by  adding  dummy 
records . 


Two  records  were  merged  if  they  had  the  same  year,  day,  hour  and 
observation  code  (column  8 or  9) , but  different  values  for  the  cells.  The 
subsequent  records  were  merged  into  the  superset  record.  If  a record  was 
repeated  more  than  once  the  extra  records  were  deleted.  In  the  case  of  miss- 
ing data,  records  with  code  9 in  column  8 or  code  3 in  column  9 were  added 
depending  upon  the  data  surrounding  the  missing  data.  (Code  9 indicates  no 
echoes  existed  within  cells  11-74.  Code  3 indicates  observation  not  available.) 
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Finally,  the  corrected  data  for  five  stations  was  put  on  a 7-track  tape. 
There  are  five  files  on  the  tape,  one  for  each  station. 


SOFTWARE  AND  ALGORITHMS 

The  program  written  is  composed  of  a main  program  and  four  sub- 
routines. In  the  main  program  parameters  are  initialized  and  a matrix  which 
contains  the  cell  values  for  8 groups  for  six  years  is  obtained.  Along  with 
this  matrix  two  vectors  are  formed;  one  for  the  days,  and  one  for  the  hours. 
The  frequency  distribution  (TFRD) , probability  distribution  (P) , normalized 
probability  distribution  (XR) , correlation  coefficients  (CC) , and  the  lag 
correlation  coefficients  (CCN)  are  calculated  for  a season.  This  process  is 
repeated  for  four  seasons.  The  program  was  run  for  all  five  stations.  The 
probability  matrix  P is  defined  by: 


66 

ZlFRD  (i,  j)  - .375 


sample  size  + .25 


K = 66,  65,  . . . , 2 


For  K = 1 (blanks) , then 


P 


(i, 


j) 


66 

ZlFRD  (i,  j)  + IFRD  (i,  j) 

i=2 2 

sample  size  + .25 


.375 


The  normalized  probability  matrix  XR  is  defined  by: 

XR(L,  j)  = C*T  - ANUM/DENOM 

where 

ANUM  = AO  + A1*T 

DENOM  = 1 + Bl*T  + B2*T**2 

AO  = 2.30753 
A1  - 0.27061 


B1  = 0.99229 
B2  - 0.04481 


C 


I 


If  P (L,  J)  ^ 0.5  C = 1 

T = SQRT  (ALOG ( 1/P (L, J ) * *2 ) ) 

If  P (L,  J)  > 0.5  C = -1 

T = SQRT  (ALOG (1/ (1-P(L,J) ) * *2 ) ) 


The  correlation  coefficient  equation  used  is: 


P = — 3 — 
xy  a a 
x y 


.hero  q - (>. 

= standard  deviation  of  group  x 
o = standard  deviation  of  group  y 


i 

I 

i 
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3. 


Conditional  Probabilities  and  Scores 


INITIATOR:  I.  GRINGORTEN 

PROBLEM  NO.:  4862B  PROJECT  NO. : 8624 


BACKGROUND 

Given  the  frequency  of  a weather  condition,  it  is  possible  to 
estimate  the  probability  of  its  occurrence  or  to  calculate  its  condi- 
tional probability  following  various  antecedents.  From  Problem  No.  4862A 
(Section  2),  a matrix  of  correlation  coefficients  was  found  between  seven 
predictors.  Similarly,  a vector  of  correlations  between  a predictand  and 
each  of  the  seven  predictors  was  formed.  The  objective  of  this  task  was 
to  use  these  correlations  to  find  the  conditional  probabilities  of  a pre- 
dictand, given  each  of  the  seven  predictors.  Also  the  ability  to  estimate 
the  conditional  probabilities  was  evaluated. 

The  purpose  of  this  task  was  (1)  to  generate  the  software  to 
calculate  the  necessary  probabilities  and  other  associated  statistics; 
and  (2)  to  calculate  the  necessary  statistics.  The  results  of  this  task 
and  the  task  described  in  Section  2 were  used  by  the  researcher  and  pub- 
lished^ . 

ANALYSIS  AND  RESULTS 

The  conditional  probability  was  computed  using  the  following  equa- 
tion : 


1 Gringorten,  I.  I.,  "Multi-Prediction  Conditional  Probabilities,"  AFGL-TR- 
76-0239. 
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P (Y  = y I X = x) 


P(Y  = y n X = x) 
P(X  = x) 


In  the  more  general  case,  if  Y is  the  predictand,  and 
are  n predictors,  then  it  is  desired  to  find 


P (Y  > Y |x, , . . . , X ) 
— cl  n 


where  Yc  is  a threshold  value  of  the  predictand. 


The  following  values,  calculated  in  Problem  4862A  (Section  2), 
were  used  in  this  analysis: 


p.  . = E[X.X.] 
ID  i D 

P.  = E [YX . ] 

i l 


where  c. . are  the  correlation  coefficients  between  predictors  X.  and  X 
ID  1 

and  o ^ are  the  correlation  coefficients  between  predictand  Y and  the  i 

predictor,  X^. 

The  assumption  of  linear  dependence  of  the  predictand  and  the  pre- 
dictors is  necessary  to  derive  the  following  analytical  model  for  conditional 
probability: 

y=a,x,+...+ax+bn  (1) 

r 1 1 n n 


ih 


where  the  a^'s  are  partial  regression  coefficients,  b is  a coefficient  with 
magnitude  such  that  the  normality  of  y is  preserved,  and  n is  random  and 
normally  distributed. 

By  squaring  both  sides  of  equation  (1)  and  calculating  expected 
values,  we  get  the  following  expression: 


9 


n 

1 = £ ai  + 2 £ 

i=l  i?<j 


By  multiplying  both  sides  of  equation  (1)  by  and  taking  expected 
values,  we  get: 


a . a . p , . + b^  (2) 

i D il 


P. 

l 


a . + 

i 


£ 

j/i 


a . p . . 

D 13 


(3) 
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By  multiplying  both  sides  of  equation  (4)  by  the  inverse  of  C,  we  have: 


A = C_1R  - CR/|cj 

where  |c|  is  the  determinant  of  matrix  C 
C is  the  adjoint  matrix  of  C 

From  equation  (2),  we  derive  the  expression: 


(5) 


or 


b = /l  - a,  p,  - . . . - a p 
11  n n 


(6) 


From  equation  (1)  for  specific  values  of  x, , . . . , x , the  value  of 

1 n 

exceed  as  frequently  as  y exceeds  an  assigned  minimum  y^ . Thus 


p(n  > nc)  = p(y  > yc|x  , 


/ x ) 
n 


(7) 


The  probability  that  n will  exceed  is  the  conditional  probability 

y will  exceed  y . The  solution  for  n is  given  by 
c c 


nc  = (yc  ' aiXl 


. - a x )/b 
n n 


(8) 


where  a^  are  given  by  equation  (5),  and  b is  given  by  equation  (6) 


n win 


that 
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Once  x and  v are  known,  then  ■ is  determined  by  eouation  (8). 
i c 

We  then  use  equation  (7)  to  compute  the  conditional  probability.  A useful 
approximation  used  in  this  analysis  is  provided  by  the  following  algorithm*: 

p(n  ) = 2 + m [2d  + c.n  + c n2  + c,n3  + c^n4)4]-1 

c Ic2c3c4c 

where 

C = 0.196854 
C2  = 0.115194 
C3  = 0.000344 

C,  = 0.019527 

4 

2 = 0,  m = 1 for  n _>  0 
2 = 1 , m = - 1 for  h < 0 

For  this  task  we  modified  the  routine  to  calculate  a matrix  or  cor- 
relation coefficients  between  seven  predictors,  and  also  a vector  of  corre- 
lations between  a predictand  and  each  of  the  seven  predictors.  The  corre- 
lations previously  calculated  were  used  to  find  the  conditional  probability 
of  a predictand,  given  each  of  the  seven  predictors.  This  predictand  is 
linearly,  but  stochastically,  dependent  on  the  predictors,  and  all  variables 
are  normalized  and  Gaussianly  distributed.  In  addition  to  the  standard 
correlation,  lag  correlations  were  calculated  between  the  seven  groups 
(groups  two  to  eight)  used  as  predictors  and  the  first  group. 

Given  the  probability  matrix  P,  and  the  correlation  coefficient 
vector  R , the  values  calculated  for  Problem  4862A  (Section  2).  addi- 

l 

tional  software  was  written  to 
1)  Calculate  P and  |p| 


* National  Bureau  of  Standard  (1964)  Handbook  of  Mathematical  Functions  with 
Formulas",  Graphs  and  Mathematical  Tables,  Applied  Mathematics  Series  No.  55 , 
U.  S.  Government  Printing  Office,  Washington,  D.  C. 


12 


fl 

2) 

Determine  A = PR/ | 

P | where 

A = 

* 

a7 

, T 

, 1/2 

3) 

Find  b = (1  - A 

R ) 

4) 

Extract  the  values 

; v ••• 

, X7  (calculated  pre- 

viously)  for  each 

hour  in 

the  sample. 

Also  extract  four  selected  values  of  Y from 

Q 

previous  results  ^nd  calculate 

Yc  - alXl  ...  - a7x? 

^c  b 

Then  calculate  the  conditional  probability 

1 


Prob  (rpnc)  = 


r 2 3 4 -i4 

2 [ i + c^n  + c2n  + c3n  + c4n  J 


for  n > o 


and 


Prob  (n>n  ) = l 
— c 


4 -.4 


2[1  + c^n  + c2n2  + c3n3  + c^n*1] 


for  n < 0 


where  = .196854 
C2  = .115194 
= .000344 


C„  = .019527 
4 


5)  Find  the  distribution  of  values  of  Prob  (n  > n ) 

— c 

6)  Determine  the  score  (C  - T)/C 

2 

where  C = E (Prob  (1,  L)  - V) 

2 

T = E (Prob  Y - V) 

and  V = 1 when  the  number  of  octas  on  the  Nth 
hour  equals  or  exceeds  the  chosen 
value  of  L 


V = 0 when  the  number  of  octas  is  less 
than  the  chosen  value  of  L 
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Prob  (1,  L)  = the  unconditional  probability 
that  the  number  of  octas  in 
Group  1 will  equal  or  exceed  L 

Prob  Y = the  conditional  probability  that 

the  number  of  octas  in  Group  1 

will  equal  or  exceed  L;  i.e., 

Prob  Y = Prob  (n  > n ) 

— c 


The  above  parameters  were  calculated  and  printed  for  four  seasons  and 
five  sites. 
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4,  In  Flight  Line  of  Sight 


INITIATOR:  E.  BERTONI 

PROBLEM  NO.:  4863  PROJECT  NO. : 8624 


BACKGROUND 

The  purpose  of  this  task  was  to  calculate  vertical  profiles  of 
the  probability  of  clear,  cloudy  and  hazy  lines  of  sight  from  aircraft, 
at  1000-foot  increments  over  selected  areas  of  the  Northern  Hemisphere. 
Seasonal  as  well  as  total  data  calculations  were  performed.  Similar  re- 
quests had  been  performed  previously  and  software  existed  from  these  tasks. 


ANALYSIS  AND  RESULTS 

We  received  four  programs  necessary  for  performing  this  task. 

However,  they  required  modifications  to  complete  the  task  as  requested 
(e.g. , to  output  the  total  number  of  observations  for  each  season).  The 
documentation  was  outdated  and  did  not  describe  the  programs  as  we  received 
them.  Hence,  a new  set  of  documentation  was  required  in  addition  to  the 
programming  changes. 

The  first  subtask  was  to  verify  the  programs.  The  programs  were 
checked  out  by  using  all  four  previous  copies  of  the  data  tapes  and  the  data 
set  on  the  private  desk  file.  The  results  were  compared  to  those  of  prior 
runs.  (Problems  occurred  and  will  be  discussed  briefly  later  in  this  section.) 
After  this  check  was  completed,  the  last  program  of  the  package  was  modified 
to  make  calculations  at  46  intervals  instead  of  the  original  eight.  This 
change  was  checked  by  summing  over  a number  of  intervals  and  comparing  the 
sum  with  the  results  of  an  eight-interval  run.  Once  this  was  done,  four 
additional  boxes  were  added  to  the  old  data,  and  the  resulting  file  was  pro- 
cessed. 
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The  new  data  was  input  to  the  first  program  (LANLAT)  which  per- 


cks 

for  valid  data, 

such  as 

1) 

1 < MONTH  < 12 

2) 

63  < YEAR  < 76 

3) 

-90  <_  LATITUDE  < 

90 

4) 

-180  <_  LONGITUDE 

< 180 

5) 

ALTITUDE  > 0 

Data  not  within  the  specified  ranges  was  rejected.  Latitude  was 
then  converted  to  0®  <_  LAT  180°  and  sorted  in  hundredths  of  a degree, 
where  0°  is  the  south  pole.  Longitude  was  converted  to  0°  LONG  < 360° 
west  longitude  and  stored  in  hundredths  of  a degree.  Ail  line  of  sight 
indicators  (+60°,  +30°,  HOR,  -30°,  -60°)  which  are  not  equal  to  1 (clear), 

2 (cloudy)  or  3 (hazy)  are  set  to  zero. 

Valid  sorted  data  was  stored  on  disk  and  then  sorted  by  longi- 
tude and  output  to  tape.  The  second  program  used  the  SORT/MERGE  utility 
program  to  merge  the  output  file  from  the  first  program  with  the  old  data 
file.  The  old  data  file  is  also  sorted  by  longitude.  (This  is  the  output 
from  the  previous  run.) 

A third  program  (CRESIS)  was  then  run,  using  the  merged  data  t< 
create  a SCOPE  INDEXED  SEQUENTIAL  (SIS)  file  with  a key  based  on  the  longi- 
tude. This  program  aborted  if  the  input  file  was  not  sorted  by  longitude. 
Since  the  operating  system  of  the  CDC  6600  had  been  changed  since  the  last 
run  by  the  addition  of  the  RECORD  MANAGER,  there  were  problems  with  gettinq 
this  program  to  run  properly.  Modifications  were  made  to  the  software  in 
order  to  run  under  the  existing  system. 

Additional  complications  arose  whi le  trying  to  process  the  data 
through  the  fourth  program  (LOSANL)  which  performed  the  actual  calculations 
for  percent  occurrence  of  cloud-free-1 ines-of-sight . Cards  were  read  to 
select  the  months,  altitudes,  and  area  of  the  world  desired  for  a particular 
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run.  The  problem  was  in  determining  what  : < -if:  input  th»  program  had 

been  designed  to  a :ej  t.  A numbei  i nodifi  it  had  t rna  ie.  For 

example , the  ngitud<  t red  ■.  r<  iegrees  ind  360  degrees , 

: ut  wr>  ib  ul.it  ion;  f-  i the  hemis: . : .•  • r>  w-  r-  r-  ; ie  • • 1,  th.  interval  had 
to  be  specified  as  0 - 359.99,  since  o and  360  wen-  treat'  d a.;  the  same  point 
In  addition  to  this,  the  : rogram  w.i  n t design-  i * t in  interval 

wl  :]  : ■ i iegrees . rhi  ua  t r , i t 

had  ■■•••’er  :r : • i , «nd  : ’ w i t • it  . • ild, 

ittempt  wet  made  t rrect  this.  Anothei  tl  prog ram 

was  that  it  could  not  accept  mor-  t nan  n<  - s<  - f da*  i a*  i * : me , despite 
inii  ations  that  it  could. 

Sub  ei t to  these  limitations,  the  program  was  run  for  20  geographic 
areas  for  four  seasons  and  also  f i th(  wh  le  tear.  printed  output  from 

these  rui  ii  lyed  th<  iltitudi  inter  il  Lnv  Lved  ind  the  number  of  obser- 
vations of  each  type  (clear,  hazy,  rloudy)  at  each  altitude  and  at  each  of 
th<  horizon  (horizontal , • i , +<  , - 30° , - • ’ ) . The  pro- 

: r eacl  f thes<  weri  printed.  pi  jram  was  modified 

sc'  the  total  number  of  observations  for  each  alt  it  id<  interval  and  overall 
intervals  were  calculated  and  displayed. 
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5.  Kwajalein  Missile  Range  Radar  Conversion 


/ 


INITIATOR:  A.  FAIRE 


PROBLEM  NO.:  4B96A  PROJECT  NO.:  6690 


BACKGROUND 

Tracking  of  missiles  is  performed  at  different  sites.  The  data  from 
these  sites  was  compared  to  obtain  the  optimum  representation  of  the  flight 
characteristics  of  the  missile.  Each  site  has  its  own  type  of  computer 
hardware  for  obtaining  the  data.  Before  the  analysis  of  the  data  could 
begin,  the  data  tape  was  converted  to  a format  compatible  with  the  computer 
that  was  to  be  used  in  the  analysis. 

Conversion  of  the  data  tapes  required  a knowledge  of  the  manner  in 
which  the  data  was  obtained.  It  was  necessary  to  be  informed  of  the  com- 
puter type,  its  word  size,  byte  size,  character  code  used  and  method  of 
storing  floating  point  arithmetic  numbers.  Also,  quite  often  in  order  to 
decrease  tape  reading  time,  data  was  blocked  with  more  than  one  record  per 
block.  The  units  of  the  data  also  had  to  be  known  before  processing  could 
begin.  Once  conversion  was  accomplished , data  was  printed  and  plotted 
for  comparison  and  further  analysis. 


ANALYSIS  AND  RESULTS 

The  radar  trajectory  data  tapes  were  processed  at  the  Honolulu 
Data  Reduction  Facility  on  an  IBM  7094  computer.  This  machine  has  36-bit 
words.  The  tapes  were  7-track  with  a density  of  800  B.P.I.  The  record 
size  was  24  words  with  a block  size  of  1 record  per  block.  Table  5-1  shows 
the  data  as  it  appeared  on  the  tape,  with  associated  units.  The  record  for- 
mat type  is  U.  Accordingly,  the  input/output  processor  of  the  computer  did 
not  search  for  record  delimiters  or  control  words;  and  if  they  did  exist, 


18 


they  would  have  been  considered  as  part  of  the  data.  In  addition,  each 
block  was  read  as  if  it  contained  one  logical  record.  The  layout  of  a re- 
cord is  shown  in  Figure  5-1.  Each  word  was  composed  of  an  exponent  (bits 
31  through  34),  its  associated  sign  (bit  35),  the  number  (bits  1 through  30) 
and  its  associated  sign  (bit  36) . This  can  be  represented  mathematically 
by  the  following: 

S E 

F = S J (2)  E (1) 

N 

where  F is  the  true  value  of  the  data  and  S , N,  S , and  E are  as  described 

N E 

in  Figure  5-1. 


In  order  to  arrive  at  the  final  converted  value,  two  important 
computing  steps  were  required.  The  first  is  to  unpack  each  36-bit  7094 
word  into  a CDC  floating  point  number.  The  CDC  floating  point  arithmetic 
is  described  in  Section  29  of  this  report. 

Processing  of  data  was  begun  by  buffering  the  24  words  (36  bits/ 
word)  into  15  CDC  60-bit  words.  It  can  readily  be  seen  that  5 IBM  36-bit 
words  equal  3 CDC  60-bit  words.  If  we  process  5 input  words  at  a time,  a 
loop  can  be  set  up  to  convert  all  the  data.  This  was  accomplished  by  shift- 
ing bits  until  each  CDC  word  contained  36  bits  of  data  (see  Figure  5-2). 

A loop  was  then  used  to  shift  the  bits  associated  with  each  component  of  the 
36-bit  word  (that  is,  number,  exponent,  sign  of  number  and  sign  of  exponent), 
separate  the  components,  and  apply  equation  (1)  to  obtain  F.  At  this  point, 
the  data  was  in  the  desired  format. 

The  SUA  CDC  6600  Radar  Tape  Format  is  shown  in  Table  5-2.  Header 
information  was  supplied  by  the  user,  and  data  was  output  to  the  printer  as 
well  as  to  tape.  Provision  was  made  for  processing  more  than  one  tape  if 
desired. 


Plots  of  the  data  were  also  obtained.  Becuase  of  the  large  volume 
of  data  being  processed,  it  was  necessary  to  plot  only  a part  of  the  data 
for  each  parameter  at  a given  time.  The  plots  gave  a representation  of  the 
raw  data  as  it  appeared  on  the  data  tape.  Plot  limits  were  set  automati- 
cally by  the  program.  The  data  was  read  once  to  determine  maximum  and  mini- 
mum values,  then  the  tape  was  rewound  for  actual  processing.  The  program 
allowed  scaling  of  the  axis  in  desired  units. 

Once  conversion  was  completed,  the  data  was  then  compatible  both 
in  format  and  units  for  use  by  all  users  at  AFGL.  Determination  of  velo- 
cities, densities,  etc.  could  then  be  determined.  Future  tapes  can  be  con- 
verted with  little,  if  any,  loss  of  time  and  a minimum  of  supervision.  It 
is  expected  that  data  to  be  processed  will  be  received  at  monthly  intervals. 

In  addition  to  the  conversion  of  this  particular  data  set,  these 
algorithms  offer  a basis  for  conversion  of  other  data  sets  produced  on 
similar  36-bit  machines. 


TAPE  DATA  RECORD  FORMAT  (U) 


1)  A negative  GMT  signifies  end  of  data  (-9999.  in  Word  2) 

2)  An  EOF  marker  signifies  end  of  tape 


TABLE  5-2  SUA  CDQ  6600  RADAR  FORMAT 


The  radar  data  type  format  consists  of  a 20  word  header  record  followed  by 
96  word  data  records.  The  96  word  data  record  contains  6 samples  of  16  words 
each.  At  the  end  of  the  header  record  and  the  96  word  data  records  is  an 
end  of  record  word. 


TABLE  5-2  SUA  CDC  6600  RADAR  FORMAT  (CONTINUED) 


16  Word  Data  Records,  6 Samples  per  96  Words  Record 

Word 

Mode 

Parameter 

Units 

1 

Floating 

Time  of  Day  (GMT) 

Milliseconds 

2 

3 

4 

5 

6 

Floating 

Slant  Range 

Yards 

7 

II 

Range  Rate 

Yds/sec 

8 

It 

SNK 

dB 

9 

II 

Azimith 

Degrees 

10 

II 

Elevation 

Degrees 

11 

12 

Integer 

Auxilliary  Word  1 

13 

Integer 

Auxilliary  Word  2 

14 

15 

16 

25 


b.  Conversion  Of  ALCOR  KMR  Metric  Data  Tape 


INITIATOR:  A.  FAIRE 

PROBLEM  NO:  4896B  PROJECT  NO:  6690 


BACKGROUND 

Four  metric  d ?ta  tapes  generated  by  Lincoln  Laboratory  in  Lexington, 
MA,  were  received  for  further  processing.  These  data  tapes  were  produced  on 
an  IBM  370  computer  but  the  data  was  formatted  to  be  compatible  with  standard 
CDC  6600  hardware.  This  data  is  similar  to  the  Kwajalein  Missile  Range 
Radar  Data  (see  Section  5) . Both  sets  of  data  describe  trajectory  information 
from  Rocket  All. 605.  The  purpose  of  this  analysis  was  to  convert  the  data  into 
the  British  Unit  System  and  to  re-format  the  data  eliminating  unused  data. 

This  data  was  made  available  for  printer  output,  and  subsequently  plotted. 
Checks  were  performed  to  determine  missing  data  and  an  asterisk  placed  as 
a flag  on  printer  output  when  this  occured. 


ANALYSIS  AND  RESULTS 

Problems  due  to  errors  in  the  documentation  arose  in  the  conversion 
of  these  data  tapes.  The  tape  was  expected  to  be  a 9-track,  800  (or  1600) 

BPI,  multi-file  tape,  with  a record  size  of  50  32-bit  words  and  blocked  50:1. 
However,  we  determined  that  the  tape  wa  ; actually  7-track,  800  BPI,  with  a 
record  size  of  50  60-bit  words  and  blocked  10:1.  Thus,  it  is  uncertain  if 
future  tapes  will  be  compatible  to  the  present  method  of  obtaining  data.  These 
differences  were  conveyed  to  Analysis  and  Simulation  Branch  personnel  who  then 
conveyed  the  information  to  Lincoln  Laboratory. 

The  contents  of  the  header  and  data  records  of  the  input  data  tape 
are  explained  in  Tables  6-1  and  6-2  respectively  as  they  apply  to  each  radar. 
The  header  record  is  almost  identical  from  radar  to  radar.  Data  records  differ 
in  content  and  definition  for  each  radar. 
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Once  the  data  could  be  read  and  converted  into  the  British  Unit 
System,  other  problems  were  encountered.  The  data  rate  was  expected  to 
be  1000  points/second,  but  was  found  to  be  200  points/second.  The  time 
changes  were  not  in  even  increments  and  occasionally  differed  by  .00003 
seconds  from  the  expected.  This  made  it  necessary  to  interpolate  the  time 
in  order  to  obtain  the  data  for  every  .1  second  as  desired  for  SUA  tape  format. 
Word  #19,  the  qualifier  word,  did  not  meet  specification.  After  discussions, 
it  was  decided  to  ignore  it.  To  find  missing  data,  the  words  were  checked 
to  see  if  they  were  zero. 

Table  5-2  in  Section  5 lists  the  format  for  SUA  CDC  Radar  Format 
used  for  this  task.  The  header  was  created  from  input  card  data.  The 
program,  regardless  of  the  number  of  files  on  the  input  tape,  generates  one 
output  tape  composed  of  all  files  of  the  input  tape. 

Plots  of  the  data  were  created  by:  (1)  plotting  180  data  points  for 
each  variable  at  one  time;  (2)  re-initiating  plot  parameters;  and  (3)  returning 
to  the  first  graph  to  plot  another  180  points.  Checks  were  made  to  insure 
that  the  data  variables  were  not  zero.  These  plots  were  used  for  comparison 
of  the  raw  data  generated  from  different  sites. 

Four  tapes  have  been  processed  according  to  the  procedure  described 
above.  Although  these  tapes  were  in  agreement,  it  may  be  necessary  to  make 
modifications  to  account  for  any  discrepancies  in  future  tapes.  This  data, 
along  with  that  received  from  KENTRON  Honolulu,  will  be  used  in  analysis  of 
density  data.  It  is  expected  that  more  tapes  will  be  recieved  for  Rocket 
#A10. 712-1  and  Rocket  #A10. 712-2  which  will  require  similar  formatting. 
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Table  6-1  Header  Record 


Def init ion 

Tape  identifier  = X 

where : 

X = 21  for  ALCOR 


File  designator  = MTXXY  ABC 
where : 

MT  = Metric  Tape 
XX  = Any  integer 

Y = Letter  identifying  the  company /computer  for  which  a 
tape  has  been  tailored. 

AM C = Three-digit  number  for  Lincoln  Laboratory  identifi- 
cation purposes. 

Mission  date  = XXYYY 


where : 

XX  = The  last  two  digits  of  the  year  of  launch 
YYY  = Julian  day 

Object  - Tracked  object  to  which  metric  data  in  this  file  applies. 
Launch  time  (GMT) 


Computer  Name  (The  original  tapes  are  generated  for  IB.'!  360/ 
370  computers.  If  a tape  is  tailored  to  a user's  computer, 
the  computer's  name  is  entered  here.) 

Program  version  date  (Each  radar  has  a separate  magnetic  tape- 
generating program  whose  version  dates  are  entered  here.) 


Table  6-1  Header  Record  (cont.) 


Word 

Def ini tion 

12 

Number  of  time  spans  (The  number  of  trackinq  intervals  or 
data  spans  included  in  this  file.  The  maximum  number  has 
been  arbitrarily  set  to  10.) 

Span  1,  start  time  (GMT'  - Beqin  time  of  first  interval. 

1 

1 14 

Span  1,  stop  time  (GMT)  - End  time  of  first  interval. 

15 

Span  1,  data  rate  (pts/sec)  = X (X  200) 

43 

Source  of  beacon/skin  separation  = X*  j 

| 

where : 

X = 1 for  the  comparison  of  beacon  and  WB  range  trajectories 

2 for  the  comparison  of  beacon  and  NB  range  trajectories 

t 

3 for  the  determination  from  a ’JB  pulse  shape 

1 

1 

4 for  the  nominal  value  in  the  LEXCAL  array  (Words  649-642) 

*The  ALCOR  data  users  manual  may  be  referenced  for  interpretation  and  units. 
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Table  6-2  Data  Record 


Word 


Definition 


1-2 


GMT  = time  at  incidence 


Time  after  launch  = TAL  at  incidence 

= GMT  - launch  time  (sec) 


where : 

GMT  = Word  1 of  MT  data  record. 

Launch  time  (sec)  = Word  8 of  MT  header  record. 


5 


6 


7 


Mode  Word  1,  range  tracking 
Status  = T (track) 

= N (not  in  track) 
Mode  = C (Centroid) 

= E ( Edge ) 

= G (Gated  edge) 
Waveform  = N (NB) 

= W (WB) 

= B (Beacon) 

Mode  Word  2,  angle  tracking 
Status  = T (track) 

= N (not  in  track) 
w 'eform  = N (NB) 

= W (WB) 

= B (Beacon) 

Range  (km) 

Azimuth  (rad.) 
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8 


Elevation  (rad.) 


Table  6-2  Data  Record  (cont.) 


Def inition 


Not  used. 


Azimuth  offset  (rad.) 


Elevation  offset  (rad.) 


Beacon/skin  separation  (km) 


13-14 


Ranqe  tropospheric  refraction  correction  (km)  = 

Elevation  tropospheric  refraction  correction  (rad.)  = Elrp 


where : 


Rt  = Ar 

eIt  = ~aei 


Not  used. 


S/N  (dB)  = Signal  (dB)  - Noise  (dB) 


Not  used. 


Qualifier  Word  - The  quality  of  each  raw  data  point  is 
checked.  Plots  of  the  raw  data  and  plots  of  the  residuals 
from  a polynominal  fit  to  the  raw  data  are  used  as  a basis. 

A bit  set  to  one  in  the  bits  listed  below  indicates  a message 
to  the  user.  If  the  message  (flag)  refers  to  a particular 
span  of  data,  it  should  not  be  used.  The  following  qualifiers 
are  reserved  for  ALCOR. 


Range 
Azimuth 
Elevation 
Range  rate 
Azimuth  offset 


Elevation  offset 


S/N  range 
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■e  6-2  Data  Record  (cent.) 

Definition 

25  Undefined 

24  Undefined 

23  Marks  the  last  data  records  of  a file.  "Zero" 

records  are  added  until  the  last  block  is  filled. 

The  end-of-file  mark  is  then  added. 

22,21  Undefined 

20  Bit  is  set  for  every  data  point  (containing  coherent 

range  rate,  its  associated  TAL,  and  intervening  zeros) 
over  the  interval  for  which  coherent  range  is 
available. 

1-19  Undefined 

2 2 1/2 

Altitude  (km)  = |r  + R + 2RR  Ell  - R 
L e e J e 

where : 

R = Range  (km)  = [Word  6 + Word  12  + Word  13  of  MT  data  record 

R = Earth  radius  (km)  = 6378.135 
e 

El  = Elevation  (rad.)  = [word  8 + Word  14  of  MT  data  record^ 

A = Azimuth  error  (rad.) 

E = Elevation  error  (rad.) 

Corrected  range  (km)  = [word  6 + Word  12  + Word  13  of  MT  data 

record^ 

The  corrected  range  is  set  to  zero  if  Word  6 is  flagged  (Bit  32 
of  Word  19  is  set) . 

The  corrected  range  is  provided  at  maximum  available  data  rate. 
Corrected  azimuth  (rad.) 

Corrected  elevation  (rad.) 


26  j GMT  = Time  at  incidence  associated  with  the  coherent 

range  rate.  By  adding  Word  2 of  the  data  record,  AXCOR's 
precision  may  be  retained. 
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Coherent  range  rate. 


7, 


Conversion  of  I'nnormalized  Composite 
Interplanetary  Plasma  Tape 


INITIATOR:  D.  SMART 
PROBLEM  NO:  4894  PROJECT  NO:  8600 


BACKGROUND 

A tape  of  solar  wind  plasma  parameters  was  received  from  NASA, 

NSS DC . This  unnormalized  composite  solar  wind  tape  contains  data  from 
numerous  space-crafts  such  as  IMP/AIMP,  Explorer,  VELA,  HEOS1 , and  0G05. 

The  data  was  generated  on  an  IBM  360  computer.  The  tape  format  is  9-track, 
1600  B.P.I.,  with  spanned  variable-blocked  binary  records.  The  logical 
record  length  is  64  bytes,  including  a 4 byte  segment  description  word  and 
fifteen  4 byte  data  words.  The  block  size  is  25604  bytes.  (Each  byte 
contains  eight  bits.) 

It  was  necessary  to  convert  the  data  on  this  tape  into  a SCOPE 
compatible  binary  tape,  with  15  words  per  record,  1 record  per  block.  Ibis 
prerequisite  was  necessary  before  any  analysis  could  be  done  on  the  data  set. 

In  order  to  accomplish  this  task,  it  was  necessary  to  have  a know- 
ledge of  the  manner  in  which  the  floating  point  values,  characters,  and 
integer  values  were  stored  in  the  IBM  360  words.  Once  this  was  established, 
conversion  was  done  by  bit  manipulation. 

ANALYSIS  AND  RESULTS 

The  format  of  one  logical  record  is  described  in  Table  7-1.  With 
some  exceptions,  the  data  is  as  submitted  by  the  experimenters,  with  no 
normalization  performed  at  NSSDC  for  mutual  consistency.  Tin-  data  is  in 
sequential  time  order,  with  as  many  records  (between  0 and  ) for  a given 
hour  as  there  were  contributing  data  sets  for  that  hour.  The  exceptions  are 
dev  ®d  below. 
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Table  7-1  SOLAR  WIND  PLASMA  LOGICAL  RECORD  FORMAT 


► 


It 

l 


DATA  ITEM 

DATA  TYPE 

1 

COMMENT 

1.  Year 

1*4 

64,  65,  66 

2.  Decimal  day  of  year 

1*4 

January  1 = day  0 

3 . Hour 

1*4 

0,  1 , 2,  ....  23 

4.  Spacecraft  ID 

1*4 

1 = HEOS  1 
3 = Vela  3 A/Ve la  3B 
5 = OGO  5 

33  = Explorer  33 

34  = Explorer  34  (IMP  F) 

35  = Explorer  35 

43  = Explorer  43  (IMP  I) 

99  = See  following  notes 

5_  Number  of  fine  time  scale  points 

1*4 

Points  in  hour  averages 

6.  Temperature 

R*4 

°K  (blanks  for  ID  = 1 and  99  record)  j 

7.  Ion  density 

R*4 

(cm)  3 (blanks  for  ID  = 99  records) 

8.  Bulk  speed 

R*4 

km/sec,  aberration  corrected 

9.  Flow  longitude 

R*4 

In  degrees;  aberration  corrected; 
positive  for  flow  from  west  of 
sun;  solar  ecliptic  coordinates; 
word  contains  blanks  for  ID  = 1, 
34,  43  and  99  records. 

10.  Flow  latitude 

R*4 

In  degrees,  positive  for  flow  from 
south  of  sun;  solar  ecliptic  coord; 
word  contains  blanks  for  ID  = 1 , 3, 
34,  43  and  99  records. 

11.  Standard  Deviation  in  temperature 

R*4 

JK;  word  contains  blanks  for  ID  = 1, 
5,  43  and  99  records. 

12.  Standard  Deviation  in  density 

R*4 

(cm)  5;  word  contains  blanks  for 
ID  - 1,  5,  43  and  99  records. 

13.  Standard  Deviation  in  Bulk  Speed 

R*4 

km/sec;  word  contains  blanks  for 
ID  = 1,  5,  43  and  99  records. 

14.  Standard  Deviation  in  Flow 
longitude 

R*4 

Degrees;  word  contains  blanks  for 
ID  = 1,  5,  34,  43  and  99  records. 

15.  Standard  Deviation  in  Flow 
latitude 

R*4 

Degrees;  word  is  non-blank  only 
for  ID  = 33  and  35. 

Vela  3 data  was  provided  as  three  hour  averages  with  their  associated 
standard  deviations  and  with  the  number  of  measurements  contributing  to  each 
3 hour  average.  These  values  have  been  assigned  to  each  of  the  three  appro- 
priate hours. 

For  Explorer  33  and  35,  temperature  and  associated  standard  devia- 
tions were  not  provided,  but  an  equivalent  parameter  (most  probable  thermal 

2 

speed  w ) and  its  standard  deviation  OW  were  provided.  The  relation  4m  w = 

o o p o 

kT  (M^  = proton  mass,  k = Boltzmann  constant)  was  used  to  obtain  temperature 

(T)  from  the  given  information.  The  relation  a = (6t/6w  ) aw  was  used  to 

T o o 

obtain  Cf  from  the  given  w and  crw  . 

T ^ o o 

To  establish  consistency  in  the  sign  conventions  used  for  the  flow 
longitude  and  latitude  angles,  (see  Table  7-1)  Explorer  33  and  35  longitude 
and  latitude  angle  values  were  multiplied  by  minus  one,  as  were  the  0G0  5 
longitude  angle  values. 

The  data  associated  with  spacecraft  ID  = 99  was  provided  by  the  Los 
Alamos  Scientific  Laboratory  solar  wind  plasma  team.  Ibe  data  comes  from 
VELAS's  2,  3,  4,  5,  6,  from  IMPS  6,  7,  8 and  from  Explorer's  43,  47,  50; 
however,  individual  data  values  are  not  linked  with  specific  spacecraft. 

The  data  supplied  to  NSSDC  consisted  of  time,  3-hour  averaged  bulk  speeds,  and 
numbers  of  find  time  scale  values  in  the  3-hour  average.  The  3-hour  averages 
have  been  distributed  over  individual  hours  (as  for  VELA  3 data  discussed  earlier.) 

The  HEOS  1 data  was  extracted  by  NSSDC  from  Diodato,  et  al., 

"Astron.  Astrophys,  Suppl.",  20,  pp.  313-362,  1975.  These  3-hour  averages 
have  been  distributed  over  individual  hours. 

PROCESSING  PROCEDURE 

The  following  is  the  procedure  used  to  process  the  data.  One 
IBM  360  word  contains  32  bits  or  4 bytes;  for  the  CDC  6600,  each  word 
contains  60  bits.  Therefore,  each  CDC  6600  word  is  equal  to  15/8  of  an  IBM 
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I 


CDC  WORD  *6 


CDC  WORD  #7 


Figure  7-1  IBM  WORD  STORAGE  ON  THE  CDC  6600 
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ZERO  FILL 


IBM  WORD  #1 


CDC  WORD  ~2 


59 


31 


0 BIT 


ZERO  FILL 


IBM  WORD  82 


CDC  WORD  83 


59 


31 


0 BIT 


ZERO  FILL 


IBM  WORD  #6402 
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0 BIT 


CDC  WORD  #6402 


Figure  7-2  RESULT  OF  UNPACKING 


360  word,  or  8 CDC  6600  words  will  represent  15  IBM  360  words.  Figure  7-1 
represents  how  these  eight  CDC  6600  words  appear  in  storage.  It  takes 
3414  CDC  6600  words  to  store  the  first  block  of  data.  An  array  of  this  size 
was  set  up  and  one  block  of  data  buffered  in.  The  subroutine  to  unpack  the 
data  was  then  called.  This  subroutine  shifts  each  32-bit  IBM  360  word  into 
one  CDC  word  of  60  bits.  An  array  was  set  up  to  accept  these  values.  The 
array  size  was  3414  x 15/8  = 6402,  rounded  to  the  nearest  integer.  This 
was  accomplished  by  setting  up  2 loops,  one  to  unpack  15  IBM  360  data  words 
at  a time,  and  one  to  increment  a counter  by  15  to  process  the  entire  block. 
Once  this  was  completed,  the  values  were  then  converted  into  compatible 
CDC  6600  integer  numbers,  as  shown  in  Figure  7-2. 

The  next  step  was  to  convert  the  IBM  360  floating  point  numbers 
from  their  integer  representation  to  CDC  6600  floating  point  representation. 
An  IBM  360  floating  point  value  was  represented  as  shown  in  Figure  7-3. 

8 0 0 7 


s 

I 

EXPONENT 

MANTISSA 

G 

N_ 

0 1 18  31  BIT 


Figure  7-3  FLOATING  POINT  WORD 


The  binary  point  of  this  word  is  just  before  bit  8,  and  the 
exponent  has  a base  of  16  associated  with  it.  It  also  has  a bias  of  64;  that 
is,  64  must  be  subtracted  from  the  exponent  before  conversion.  Thus,  mathe- 
matically the  number  may  be  represented  as 


x 


(S)  (M)  (16) 


EX-64 


where 


S = Sign,  M = Mantissa,  and  EX  = Exponent 
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For  a description  of  the  CDC  6600  floating  point  arithmetic,  see 
section  29. 


To  convert  each  IBM  360  floating  point  number,  the  bits  of  each 

word  segment  were  shifted  into  separate  words.  Since  the  binary  point  is 

-24 

to  the  left  of  the  mantissa,  the  mantissa  was  multiplied  by  2 to  shift 
the  decimal  point  to  the  right.  The  formula  for  conversion  to  a CDC  6600 
floating  point  word  is 


x 


(2_24)  SM  (16) 


EX-64 


In  this  data  set,  words  7 through  16  of  each  record  were  converted 
to  floating  point.  This  was  accomplished  by  setting  up  a loop  to  call  a 
subroutine  which  does  the  conversion  as  described  above.  After  completion  of 
the  loop  a record  was  written,  and  the  next  record  of  the  block  was  processed. 
A new  block  of  data  was  buffered  in  after  the  400  records  of  the  block  were 
processed  and  written. 

When  an  end-of-file  was  reached,  a check  was  made  to  detemine  if 
all  files  were  processed.  If  processing  was  to  continue,  the  program 
returned  to  buffer  in  more  data.  If  processing  was  to  stop,  a message  was 
printed  indicating  the  last  file  processed,  and  the  program  then  terminated. 

The  general  method  used  to  convert  this  tape  may  be  cm;  loyed  in 

the  future,  with  modifications  made  to  block  length  and  size  of  record  to 
correct  future  tapes.  This  tape  was  submitted  for  approval  and  accepted.  No 
further  analysis  was  done  at  that  time. 


40 


3. 


Conversion  Of  OPAQUE  Data  Tapes 


INITIATOR:  E.  SHETTLE 


PROBLEM  NO:  4907  PROJECT  NO:  7621 


BACKGROUND 

The  purpose  of  this  task  was  to  generate  software  to  decode  data 
received  from  NATO's  Atmospheric  Optics  Measurement  Program  (OPAQUE) . This 
data  is  in  American  Standard  Code  for  Information  Interchange (ASCII ) 
character  format  on  paper  tape.  A sample  paper  tape  was  supplied.  The 
paper  tape  was  decoded  and  the  data  put  on  a magnetic  tape  with  each  block 
containing  three  logical  records.  This  magnetic  tape  was  then  used  to 
transform  the  packed  data  values  into  their  true  values  by  one  of  the  con- 
version factors  and  codes  provided.  After  the  true  values  were  successfully 
obtained,  these  values  were  then  converted  back  to  their  original  form  and 
placed  on  a second  tape.  A program  was  written  to  read  the  tapes.  The 
first  tape  generated  was  sent  to  the  OPAQUE  data  bank  in  Germany  for  veri- 
fication. The  software  was  made  flexible  so  that  if  changes  in  blocking  or 
record  size  were  made,  this  could  easily  be  incorporated  into  the  program. 

ANALYSIS  AND  RESULTS 

The  data  was  punched  on  8-channel  paper  tape  in  American  Standard 
Code  for  Information  Interchange  (ASCII) . Each  channel  transfers  into  one 
bit,  and  the  8 channels  represent  one  frame.  Each  frame  was  converted  on 
input  from  ASCII  to  internal  BCD.  The  data  was  buffered-in,  IS  words  (ISO 
characters)  at  a time.  A carriage  return  character  and  new  line  character 
were  the  only  separations  between  data. 

Figure  8-1  is  a representation  of  one  block  of  data  as  it  appears 
on  the  paper  tape.  This  data  was  represented  in  core  memory  as  one  continuous 
record.  Each  word  buffered-in  was  decoded  into  ten  words  containing  one 

v a r and  each  was  compared  with  a carriage  return  and  new  line  charac- 
t >o'.  cb  ’cter  was  saved  until  a carriage  return  and  new  line  charac- 

ter W!  r ached.  The  saved  words  represented  one  line  of  data.  These  words 
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were  decoded  according  to  the  format  in  Figure  8-2  (required  for  OPAQUE) 
with  each  output  word  representing  a data  item.  After  decoding  one  line, 
comparison  of  characters  continued  until  either  a new  set  of  carriage 
returns  and  line  feed  characters  were  reached,  or  until  the  15  words  were 
processed.  When  the  15  words  were  processed,  15  more  words  were  buffered-in, 
and  the  process  continued.  Output  consisted  of  63  data  words  and  was  blocked 
from  one  to  three  logical  records  per  block.  Output  was  desired  on  both 
printer  and  tape. 

Once  the  sample  tape  was  generated,  a second  program  was  developed 
to  read  the  generated  tape  and  print  out  the  data  with  headings  for  each  data 
item  for  identification.  The  program  simply  read  the  binary  words  into  an 
array,  and  set  up  a loop  to  output  the  words  to  the  printer  along  with  headings, 
one  record  at  a time. 

Another  program  was  developed  which  converted  the  scaled  values 

obtained  from  the  tape  into  their  true  values  by  means  of  conversion  codes 

set  up  when  the  OPAQUE  format  was  established.  Each  data  item,  as  can  be 

seen  from  Figure  8-2  (not  including  the  meteorological  parameter,  identified 

in  the  records  by  M,  which  begins  the  data  record)  is  composed  of  five 

characters.  A decimal  point  follows  the  third  character;  the  fourth  is 

the  exponent;  and  the  fifth  is  the  reliability.  Negative  exponents  were 

avoided  by  use  of  factors  which  could  be  added  to  the  exponent.  For  example, 

4 + D 

a number  12345  is  represented  as  123.0X10  where  5 is  the  reliability  and 
D is  as  follows: 


Extinction  coefficients  D = -6 
Horizontal  illuminance  D « -6 
Path  luminance  at  night  D * -8 
Path  luminance  at  day  D - -8 
Eppley  filter  order  D = -6 
Banes  filter  code  D = -6 
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The  meteorological  parameters  are  coded  as  follows: 

Word  #59  M is  used  to  indicate  start  of  meteorological  data 

Word  #59  9 - cloud  cover  value  in  eighths  with  9 used 

to  indicate  no  measurement. 

Word  #59  (32)  - wind  direction  at  10  m,  ^he  integer 

value  indicating  the  nearest  10  value  from 
true  North. 

8 = Situation  outside  OPAQUE  scheme 

9 = No  observation 

13  = Rainfall 

mm/hr  in  the  past  hour 
No  observation  or  failure:  999 

A loop  was  set  up  in  the  program  to  process  each  of  the  three 
records  in  the  block.  The  filter  options  were  read  in  from  cards  and  listed 
on  the  printer.  The  converted  data  was  written  on  a magnetic  tape  in  binary 
and  listed  on  the  printer  with  headings. 

The  last  of  the  group  of  programs  reads  the  tape  which  contained 
the  true  values,  and  converted  it  back  to  the  original  packed  data  words 
by  taking  out  the  conversion  factors  and  codes,  and  encoding  it  to  an  output 
array.  The  format  of  this  output  is  identical  to  the  first  portion  of  this 
analysis.  This  program  was  necessary  in  order  to  convert  data  obtained  from 
the  AFGL/OPA  group  into  the  tape  format  compatible  with  the  OPAQUE  Data  Bank 
This  data  could  then  be  read  by  OPAQUE. 

SUMMARY 

Four  different  programs  comprise  this'  package  of  software.  These 
programs  make  it  possible  to  convert  the  data  from  paper  tape  to  a magnetic 
tape,  from  packed  values  to  actual  values,  and  from  actual  values  to  packed 
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data.  Also,  generated  tapes  may  be  read  and  printed.  It  is  expected  that 
some  data  will  be  generated  at  AFGL  which  must  be  converted  into  packed 
values.  These  tapes  will  be  sent  to  OPAQUE  Data  Bank  to  be  concatenated 
with  data  from  other  countries  and  a new  tape  generated  from  this  new  data 
file.  Also  tapes  will  be  received  for  conversion  into  the  actual  values 
from  the  packed  format. 
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cc 

C24 1 
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4168141681416814168100911481204812458115481 
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70651706517065170651004889516595115831 

00000000000000000000000000000000000000000000000000 

0000000000000000000000000 

00000 

M9320200000201 927180000000 
R140007 


Figure  8-1  PARTIAL  LISTING  OF  PAPER  TAPE 
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WORD 

WORD 

WORD 

WORD 

WORD 
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WORD 
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WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 
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WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 

WORD 


WORD 

WORD 

WORD 

WORD 


§1  STATION  No.  (171)  Not  on  Paper  I -apt 

• 2 DATE,  Year,  Month,  Day  i 76  'I*  i N<  * n Paper  Tape 

to  be  added 

*3  TI ME  (1  300)  Not  on  paper  tap*  r b*  t ided  ir  increments  of 

100.  Ranqe  (0000  to  i ) 

• 4 DURATION  OF  MEASUREMENT  YCLI  ' 7)  Not  • »ri  Fir  t Record 

■ 5 THRU  WORD  H9  COMMENTS  IF  (CC)  N-  *mm*-nts 

•10  SCATTERING  X 100  ♦ FILTER  OPTI  >N  X 1 ♦ HUMIDITY  (D24i> 

• 11  SCATTER  TNi.  BEG  I NN I N< . Not  >n  thi  I ap.-r  ?af  « , bur  must  be 

allowed  for  (00000) 

•12  SCATTERING  FINAL  Not  on  this  Pap<  r Tap-,  but  mu.'  be  allowed 
for  (00000) 

■ 13  SCATTERING  MAXIMUM  Not  on  thi  Paper  Tap  , t>u'  must  be 

allowed  for  (OQf)OD) 

• 14  SCATTER  IN  . MINIMUM  Not  on  th.  . Paper  Tape,  but  mu  .r  be 

allowed  for  < < 0 

■ IS  NUMBER  OF  VALUES  Not  * i » ip  t : ..  , t u?  mu  • r • allow*  i f r 

• 16  EXTINCTION  ‘ill  : : v:  !*!  ie. 

•17  EXTINCTION  COEFFICIENT  (INAL  (86421) 


• IB  EXTINCTION  EFF ! I ENT  MAXIM'”  (8».42  1) 

•19  EXTI NCTIOM  COEFFICIENT  MINIMI”  (86121) 

•20  NUMBER  OF  VALUES  ( 8) 

•21  HORIZONTAL  ILLUMINANCE  BE  .INNING  (41681) 
•22  HORIZONTAL  ILLUMINANCE  FINAL  (41681) 

•23  HORIZONTAL  ILLUMINANCE  MAXIMUM  (41681) 

•24  HORIZONTAL  ILLUMINANCE  MINIMUM  41681) 

■25  NUMBER  OF  VAI "ES  (009) 
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•40  PATH  LUMINANCE  AT  DAY  (SOUTH)  (88951) 

•41  PATH  LUMINANCE  AT  DAY  (WEST)  (65951) 

•42  PATH  LUMINANCE  AT  DAY  (NORTH)  (15831) 
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•6)  RAINRATE  (000),  GROUND  CODE  (0) 


T 


Figure  8-2  REQUIRED  OUTPUT  FORMAT 
46 


d 


9.  Conversion  Of  QPAQIIF  Data  Tapes  (Modification) 


INITIATOR:  E.  SHETTLE 

PROBLEM  NO:  4911  PROJECT  NO:  7621 


BACKGROUND 

The  purpose  of  this  task  was  to  modify  the  software  generated  for 
problem  4907  (Section  8) , to  accept  slightly  different  input  data  and  to 
process  the  data. 

ANALYSIS  AND  RESULTS 

Additional  paper  tapes  were  received  that  contained  blanks  between 
each  data  item  as  shown  in  Figure  9-1.  It  was  necessary  to  create  tapes 
similar  to  those  previously  described  (Section  8 of  this  report) , but  with 
the  omission  of  the  blanks.  This  was  done  by  skipping  a character  whenever 
the  blank  was  expected  after  processing  the  data  item. 

In  addition,  it  was  desired  that  these  data  items  be  output  to 
the  tape  in  packed  format.  Hence,  instead  of  having  one  word  for  each 
data  item,  370  characters  comprised  one  record.  Therefore,  there  are 
37  words/record  and  111  words  per  block.  This  tape  was  sent  to  the 
OPAQUE  data  bank  for  verification  of  acceptability. 
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Figure  9-1  Partial  Listing  of  New  Paper  Tape 


I 

I 


1 

i 

I 


48 


4 


t 

10.  Breakdown  Problems  Associated  With  Satellites 


INITIATOR:  P.  ROTOWELL 

PROBLEM  NO:  4903  PROJECT  NO:  7661 


BACKGROUND 

In  support  of  analyses  to  investigate  the  electomagnet ic  problems 
associated  with  magnetospher ic  fluxes  incident  on  synchronoous  satellites, 
ASBC  condu  d a literature  search.  The  objectives  were  to  locate,  organ- 
ize anu  evaluate  reference  material  relevant  to: 


The  breakdown  and/or  excitation  of  plasma  external 
to  the  satellite;  and 

The  breakdown  of  insulator-coated  conducting  surfaces 
subjected  to  the  radiation  of  the  magnetospheric  environ- 
ment . 

ANALYSIS  AND  RESULTS 

The  literature  search  resulted  in  two  reports.  The  first  is 
"Electronic  Radiation  in  the  Vicinity  of  Synchronous  Orbit  Satellites: 
Literature  Search",  AFGL-TR-77-0024  . The  contents  are  briefly  described 
in  the  report's  abstract  as  follows: 

I 

I 

A categorized  list  of  representative  references  is 
presented  that  applies  to  the  task  of  estimating  and  pre- 
dicting undesireable  radiation  incident  to  a satellite 
in  synchronous  orbit.  The  categories  on  references  ad- 
dress basic  elements  of  the  task.  That  is,  some  references 
apply  to  estimating  the  particle  and  field  environment 
about  the  satellite  and  other  references  are  relevant  to 
estimating  the  behavior  of  the  medium,  either  as  support  for 

electrical  discharge  or  the  excitation  of  radiation 
through  instability  and/or  the  presence  of  plasma  bound- 
aries. Both  theoretical  and  experimental  references  are 
included  in  each  category  where  possible. 
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The  breakdown  of  plasma,  which  is  a "necessary"  prere- 
quisite to  the  generation  of  discharge  currents  in  plasma, 
Is  a much  studied  phenomenon  both  theoretically  and  ex- 
perimentally. The  plasma  breakdown  category  of  referen- 
ces, together  with  the  bibliographies  associated  with  each 
reference,  provide  a representative  information  base  to 
support  a critical  review  of  important  elements  of 
breakdown  and  electrical  discharge  literature  from  about 
1960  to  the  present  time. 


The  second  report,  concerned  with  insulator  breakdown,  is 
"Breakdown  of  Insulators:  Literature  Search"  (To  be  published) . 
The  following  abstract  briefly  summarizes  this  report: 

The  irradiation  of  dielectric  coated  metal  surfaces  by 
incident  magnetospheric  fluxes  generates  problems  related 
to  the  charging,  breakdown  and  discharging  of  the  dielec- 
tric insulator.  The  overall  search  objective  is  to  locate 
and  organize  reference  materials  that  address  significant 
aspects  of  the  problem.  The  study  divides  into  ten  some- 
what independent  subjects . 

Emphasis  is  placed  on  the  readily  available  literature 
during  the  period  up  to  and  including  1961.  The  experi- 
mental measurements  during  this  period  were  conducted  under 
idealized  conditions  to  achieve  reproductibility . That  is, 
the  breakdown  E field  is  externally  generated  and  surface  break- 
down and  premature  breakdown  through  voids  are  suppressed. 

To  unify  material,  certain  basic  references  concerned 
with  solid  state,  transport  and  kinetic  theory  are  cited 
where  appropriate.  A few  references  investigating  recent 
breakdown  theory  are  cited  and  briefly  discussed.  The 
bulk  of  the  theoretical  breakdown  references  and  all  of 
the  experimental  references  were  cited  for  the  time  period 
prior  to  1961. 


PRELIMINARY  ANALYSIS 


A preliminary  analysis  using  an  idealized  plasma  model,  break- 
down criteria  and  charge-discharge  model  was  mathematically  formulated, 
but  not  implemented. 
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Reduction  Of  Meppen,  FRG  Meteorolog ical  Data 


11. 


INITIATOR:  E.  SHETTLE 

PROBLEM  NO:  4866  PROJECT  NO:  7621 


BACKGROUND 

Meteorological  data  was  obtained  from  a field  station  in  Northern 
Germany  in  support  of  the  atmospheric  optics  measurements  being  conducted  by 
AFGL/OPA.  The  data  was  supplied  on  magnetic  tape.  The  description  of  the 
tape  is  as  follows.  The  tapes  are  9-track  and  use  IBM  compatible  format.  The 
meteorological  data  is  packed  into  20  words  per  data  set  with  8 characters/word, 
six  bits/character,  and  108  data  sets  per  block. 

The  purpose  of  this  task  was  to  read  and  decode  the  meteorological 
data  tapes,  and  to  obtain  a complete  listing  of  the  meteorological  data. 

Then,  for  particular  times  supplied  as  input,  specific  meteorological  data 
was  selected  from  one  of  the  towers.  Data  from  two  meteorological  towers 
were  available  on  the  tapes.  This  data  was  made  available  to  a calling  pro- 
gram, and  was  output  to  the  printer,  punch,  or  tape.  The  data  includes:  the 
temperature  and  dewpoint  at  2 meters;  wind  direction;  and  wind  speed,  inter- 
polated at  a height  of  11  meters. 

ANALYSIS  AND  RESULTS 

The  following  is  the  list  of  information  available  on  the  data  file. 
Information  was  updated  each  ten  minutes. 

- Identification 

- Year,  month,  day,  hour,  minutes,  change  of  measuring  procedure 

- Data  from  Tower  A 

Wind  direction  at  16m  height 
Dewpoint  at  2m  and  80m  heights 

Wind  speed  at  1 , 2,  4,  8,  16,  32,  64,  85m  height 
Temperatures  at  1 , 2,  16,  48,  80m' heights 

- Data  from  Tower  B 

Wind  speeds  at  1,  2,  4,  8,  16,  32,  64,  80m  heights 
Temperatures  at  2,  16,  48,  80m  heights 
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The  program  buffered  in  one  block  of  data  (17,280  characters)  from 
the  tape  and  decoded  it  onto  a 7-track  tape,  80  characters  per  record.  A check 
was  made  for  an  end  of  file  at  the  beginning  and  end  of  the  tape.  As  many 
files  as  are  on  the  tape  can  be  created  on  the  new  tape.  This  process  served 
two  purposes:  one,  it  created  a copy  of  the  original  tape  for  backup:  two, 
the  records  were  in  a form  which  can  be  processed  more  easily. 

Since  data  on  the  tape  was  for  both  towers,  disk  files  were  used  to 
separate  the  data  related  to  each  tower.  These  files  were  later  rewound  and 
copied  to  the  output  file  for  printing.  Headers  were  first  written  onto  each 
file.  Then  the  first  2 records  of  the  input  tape  were  read  and  a counter  was 
incremented.  The  data  was  then  written  out  onto  each  appropriate  tape.  This 
was  continued  until  50  records  were  read  and  written  for  each  tower.  The  counter 
was  reset  to  zero,  and  a new  header  was  written  on  the  files  along  with  a 
carriage  control  to  cause  the  printing  to  begin  at  the  top  of  the  next  page. 

A page  counter  was  also  inserted  for  ease  of  handling.  When  an  end  of  file 
was  reached,  the  process  terminated.  The  disk  files  were  then  rewound,  and 
copied  out  to  the  printer. 

To  obtain  specific  meteorological  data  for  a desired  time  span,  a 
new  computer  program  was  developed  which  would  read  from  an  input  card  the 
tower  (A  or  B) , starting  date  and  time,  and  ending  date  and  time  for  the 
data  desired.  The  input  parameters  were  first  printed  to  assure  the  user  of 
the  dates  he  chose.  The  input  tape  was  read  one  block  at  a time,  and  a loop 
was  set  up  to  search  for  the  starting  date  and  time.  When  it  was  present, 
the  wind  speed  at  a 10  meter  height  was  linearly  interpolated  and  the  rest 
of  the  desired  data  was  placed  in  an  array  to  be  printed  or  punched  on  cards, 
as  selected  by  the  user  as  an  option.  This  process  continued  until  the  ending 
date  and  time  was  reached,  whereupon  the  program  terminated.  The  data  was  in- 
clusive, so  the  meteorological  parameter  for  the  ending  date  and  time  was 
included  in  the  output.  To  obtain  data  for  one  specific  date  and  time,  the 
starting  and  ending  dates  must  be  the  same. 

If  data  is  desired  for  more  than  one  range  of  time  points,  new  input 
cards  must  be  read.  They  do  not  have  to  be  in  sequential  date  and  time  order. 
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These  programs  were  run  and  processed  test  tapes . These  pro- 
grams are  now  avaiable  to  process  future  meteorological  data  tapes  as 
they  are  received. 


12,  Modification  Of 
Satellite  Ephemeris  Program  (Part  I) 


INITIATOR:  J.  KLOBUCHAR 


PROBLEM  NO:  4887  PROJECT  NO:  4643 


BACKGROUND 

The  Trans-ionospher ic  Propagation  Branch  works  with  the  Air 
Weather  Service  Global  Weather  Central  (AWS/AFGWC)  in  providing  real 
time  Total  Electron  Content  (TEC)  measurements  for  military  systems  users. 
The  satellites  used  for  these  real  time  measurements  are  in  near  geosta- 
tionary orbits.  Because  of  the  remote  locations  of  the  observing  sites, 
computer  sheets  of  the  conversion  tables  from  measured  polarization  values 
to  final  TEC  values  must  be  prepared  up  to  four  or  five  months  in  advance. 
The  problem  presented  to  ASEC  was  to  provide  AWS  with  an  accurate  program 
compatible  with  their  UN'IVAC  1108. 

ANALYSIS  AND  RESULTS 

The  task  was  completed  by  merging  two  existing  programs,  namely 
ROPLOK  and  L0KANG2.  (LGKANG2  will  be  described  in  this  section.  ROPLOK 
will  be  discussed  in  Section  13.)  The  input/output  (I/O)  statements  and 
logical  steps  of  these  programs  were  modified  to  meet  the  requirements 
of  the  UN I VAC  1108  to  be  used. 

Initially  this  task  was  to  be  performed  by  generating  the 
program  LOKANG2.  LOKANG2  is  a combination  of  a satellite  ephemeris  pro- 
gram called  LOKANGL  and  the  NASA  program  MGEOS . The  resultant  program  cal- 
culates a table  of  M-values,  look  angles  for  a geostationary  satellite  as 
a function  of  satellite  latitude  and  longitude,  and  a Faraday- to-TEC 
conversion  table.  These  tables  can  be  used  by  AWS  observers  as  hourly  look- 
up tables.  While  this  program  was  not  the  final  version  sent  to  AWS,  many 
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of  the  the  techniques  used,  as  well  as  the  M-value  generation  portion, 
were  used  in  the  final  version. 

The  program  LOKANG2  was  adapted  to  accept  only  a specific 
5-card  input  deck.  The  original  program  accepted  other  inputs  as  well. 
Program  LOKANG2  was  optimized  with  the  5-card  input  in  mind.  All  un- 
necessary subroutines  were  eliminated.  Some  were  considerably  shortened 
and  included  in  their  calling  routines.  The  entire  program  was  thoroughly 
combed  for  duplication  of  efforts. 

The  purpose  of  this  version  of  the  program  is  to  enable  use  of 
L0KANG2  on  the  UNIVAC  1108  at  the  AWS-AFGWC.  To  do  this,  all  input/output 
is  done  by  calling  GREAD  and  GPRINT  subroutines  which  are  compatible  on 
UNIVAC  1108.  For  a test  case,  these  subroutines  are  simulated  on  CDC  6600 
using  Buffer  In  and  Buffer  Out  statements.* 

The  resultant  proqram  consists  of  the  following  main  section.  The 
description  of  the  routines  gives  a functional  description  of  the  program. 


NAME 


I - EPHEMERIS  SECTION 
DESCRIPTION 


LOKANG2 

WRSTA 

MNTOSC 

SHPRDC 

KEPECN 

SPSETC 


Main  program  - reads  input  data  - calls 
necessary  routines. 

Reads  M-Value  and  TEC  Conversion  Table 
output  options  and  prints  predictions. 
Extracts  predictions  from  binary  tape. 

Converts  mean  elements  to  position- 
velocity  . 

Computes  short  periodic  terms. 

Converts  mean  anomaly  to  eccentric  anomaly. 

Sets  constants  and  program  initializations 
for  all  routines. 


*For  further  description  of  LOKANG2,  see  Johanson,  J.,  "Satellite 
Beacon  Studies,"  Emmanuel  College,  Physics  Research  Division,  Nov.  1976. 
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NAME 

SPRFCO 

NEW 

CRMXA 

SPTRAN 

SPROU 

MULTIP 

ANGLE1 

ANGLE2 


I - EPHEMERIS  SECTION  (Continued) 
DESCRIPTION 

Computes  refraction  correction  due  to 
atmosphere . 

Converts  units  of  output  data. 

Prints  heading  on  each  output  page. 

Transformation  matrix  and  station  coordinates 
routine. 

Prints  out  BCD. 

Performs  matrix  multiplication 

Reduces  an  angle  to  interval  -it  to  +tt. 

Reduces  an  angle  to  interval  0 to  2n. 


II  - M FACTOR  (MAGNETIC  FIELD)  SECTION 
NAME  DESCRIPTION 

MGEOS  Calculates  M-Value  as  function  of  satellite 

positions . 

GEO  Sets  initializations  and  performs  qeometric 

calculations . 


FIELDG 

FIELD 


Calculates  geomagnetic  field,  updates  the 
POGO  (geomagnetic  field)  coefficients  to 
time  of  observation  and  prints  out  the  M-Value 
table . 


Ill  - TEC  CONVERSION  TABLE  ROUTINE 


NAME 

DESCRIPTION 

TABLE 

Produces  the  following 

output : 

1) 

TEC  (Constant  M)  to 
cor  version  table 

TEC 

(Varyinq  M) 

2) 

Faraday  rotation  to 

TEC 

conversion 

table 

In  order  to  obtain  the  above  program,  the  original  software 
was  studied  and  the  following  changes  made. 
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OPTIMIZATION  AND  CHANGES 


SUBROUTINE 

ACTION 

REASON 

CONVRT 

Deleted 

This  subroutine  was  called  only  from 
the  main  program  if  a 5-card  input 
was  not  used.  (IXYZ  / 5) 

DELEM 

Deleted 

The  only  paths  to  this  subroutine  from 
the  main  program  are  taken  if 
a one  (l)in  the  input  is  used.  No  other 
calls  are  made  to  DELEM. 

DREV 

Deleted 

DREV  is  called  only  by  DELEM  (deleted 
above) . 

NUMINT 

Deleted 

NUMINT  is  called  only  by  DREV  (deleted 
above) . 

INTGND 

Deleted 

INTGND  is  called  only  by  NUMINT 
(deleted  above) . 

DENSEL 

Deleted 

DENSEL  is  called  only  by  DELEM  and 
INTGND,  which  are  both  deleted  above. 

OSCTMN 

Deleted 

OSCTMN  is  called  only  by  DELEM  (deleted 
above) . 

PVTELM 

Deleted 

PVTELM  is  called  only  by  OSCTMN 
(deleted  above) . 

MNTOSC 

ELMTPV 

} 

Combined 

ELMTPV  is  called  only  by  MNTOSC  and 
essentially  does  the  calculations 
supposedly  done  by  MNTOSC. 

ANGLE1 

} 

Combined  via 

These  Function  subprograms  are  very 

ANGLE2 

an  Entry  State- 

similar and  were  therefore  combined 

ment 

into  one  subprogram. 

SPTRAN 

Shortened 

Several  lines  were  deleted  since  there 
was  no  path  to  them. 

SPROU 

Shortened 

With  the  5-card  input,  certain  "variables" 

in  LOKANGL  become  constants.  (Example: 
IDV  = IPOS  = IELE  = ISTA  = 0)  Therefore, 
all  IF  statements  (and  their  references) 
which  involve  these  variables  were 
deleted.  Subroutine  SPROU  also  prints 
a table  which  was  not  discussed  in 
LOKANG2  (per  discussion  with  J.  Johanson 
8/10/76) . Therefore,  the  lines  which 
generate  the  table  were  deleted. 
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OPTIMIZATION  AND  CHANGES 

(Continued) 

SUBROUTINE 

ACTION 

REASON 

JDATE 

Deleted 

JDATE  was  referenced  several  times  with 
constant  parameters.  Therefore,  it  was 
included  in  each  calling  program  (usually 
as  one  line)  instead  of  the  line  CALL 
JDATE  ( ) and  JDATt  itself  was  deleted 

CDATE 

Deleted 

CDATE  was  referenced  only  by  SPROU 

and  therefore  was  incorporated  into  SPROU. 

WRSTA  • 
WRTAP 

Combined 

via  These  two  subroutines  are  essentially 

Entry  Statement  the  same.  Since  many  of  the  calculations 

were  duplicated  in  both  subroutines, 
they  were  combined. 

TABLE 

Shortened 

Approximately  50  lines  were  deleted 
since  there  was  no  path  to  them. 

MGEOS  , 
GEOMGI 

Combined 

Since  MGEOS  was  the  only  subroutine 
to  call  GEOMGI,  GEOMGI  was  incorporated 
into  MGEOS . 

SINDF  , 

Deleted 

Since  subroutine  GEO  was  the  only  place 

COSDF  * 

in  which  these  subfunctions  were  used, 
they  were  placed  at  the  beginning  of 
GEO  and  defined  as  functions  there. 

HEAD  , 

HEADER 

Combined 

Subroutine  HEADER  is  called  only  by 
HEAD,  and  therefore,  has  been  included 
in  HEAD. 

ARSIN  , 
ARCOS  ' 

Deleted 

These  subprogram  functions  are  acceptable 
names  on  the  CDC  6600  and  are  therefore 
deleted . 

The  following  subroutines  are 

essentially  unchanged: 

SHPRDC 

NEW 

KEPECN 

CRMXA 

SPSECT 

MULTI P 

QUAD 

GEO 

AATAN 

XYZ 

SPRFCO 

AZIM 

ARCCHK 

FIELDG 

HOUk 

FIELD 
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All  subroutines  were  modified  where  appropriate  to  delete  any  un- 
used variables  or  statement  labels.  All  unnecessary  statements  were  deleted. 
The  COMMON  block  was  modified  to  reflect  usage  of  the  variables.  Certain 
programming  errors  were  corrected.  For  example:  in  subroutine  AATAN , a 

card  was  changed  to  read  1 RESULT  = 1.5707963.  In  subroutine  SPROU,  AXRATE 
was  punched  instead  of  AZRATE  and  DERATE  instead  of  DCRATE . These  were 
corrected . 

Array  DVP  was  deleted  and  replaced  by  dimensioning  DV(6)  instead 
of  DV(3).  In  the  original  program,  DV(3)  and  DVP(3)  are  placed  in  COMMON, 
causing  all  references  to  DV(4),  DV(5),  or  DV(6)  to  pick  up  the  values  in 
DVP ( 1 ) , DVP (2)  or  DVP  (3). 

In  LCKANG2  all  the  information  which  was  written  on  tapes  in 
LOKANGL  has  been  put  into  dummy  arrays  (e.g.  XTAP2 , XTAP3,  XTAP4 ) . These 
arrays  have  been  put  into  common  statements  and  inserted  into  the  proper 
routines . 


A new  variable  (HEIGT)  has  been  read  into  L0KANG2  which  is  a con- 
stant in  LOKANGL.  This  new  variable  is  the  M factor  which  was  set  equal  to 
a constant  in  LOKANGL.  This  variable  is  read  in  subroutine  WRSTA,  and 
printed  out  in  subroutine  HEAD. 

Some  of  the  output  which  was  first  printed  on  tapes  then  rewound 
and  printed  out  in  LOKANGL  has  been  stored  in  dummy  arrays  (XTP1,  XTP5)  and 
then  printed  out  in  subroutine  WRSTA. 

Further  corrections  were  needed  and  changes  were  incorporated 
after  conversations  with  AFGL/PHP  personnel.  In  subroutine  HOUR  one  of 
the  statements  has  been  changed.  The  original  statement  ZX  = 9 has  been 
changed  to  ZX  = Z.  In  the  output  page  of  LOKANG2  a new  column  has  been 
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added,  TEC/PI.  In  LOKANGL,  one  of  the  output  pages  has  two  iderfical 
columns  (RTASC  = DFCL) . The  second  column  has  been  changed  to  its  correct 
values.  Also,  in  the  output  section,  two  more  columns  have  been  changed. 

In  the  chart,  "Prepared  by  For  the  altered  column  is  E.  Longitude. 

The  original  program  was  W.  Longitude.  E.  Longitude  = 360  degrees  - 
W.  Longitude.  In  the  chart  "Look  Angles  for  ...,"  the  value  under  E.  Long- 
itude has  been  changed.  E.  Longitude  (new)  = 360  degrees  + E.  Longitude  (old) 


1108  SIMULATION 


In  order  to  run  Program  LOKANG2  on  the  1108  machine,  certain  other 
changes  were  made  to  LOKANGL.  For  example,  the  POGO  coefficients,  which 
were  read  from  permanent  file  in  LOKANGL,  are  now  read  from  a data  deck  in- 
cluded at  the  end  of  LOKANG2 . This  allows  the  program  to  be  run  indepen- 
dently of  permanent  files  located  at  only  one  facility. 

Also,  all  I/O  is  done  via  BUFFER  IN  (GREAD)  and  BUFFER  OUT  (GPRINT) 
statemei.ts  so  that  the  format  is  compatible  with  the  I/O  subroutine  which 
i utilized  at  AWS-AFGWC.  Finally,  the  card  deck  was  punched  by  an  026  key- 
punch, again  to  be  compatible  with  AWS-AFGWC. 
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13.  Modification  of 
Satellite  Fphemeris  Program  (Part  II) 


INITIATOR:  J.  KLOBUCHAR 


PROBLEM  NO:  4912  PROJECT  NO:  4643 


BACKGROUND 

Before  the  program  LOKANG2  was  sent  to  AWS , it  was  determined  that 
the  original  LOKANGL  program  had  problems  in  its  ability  to  accurately 
predict  the  satellite  longitudinal  drift.  The  problem  was  traced  to  high 
order  terms  in  the  earth's  gravitational  field.  LOKANGL  did  not  predict 
correctly  the  changes  in  satellite  position  over  extended  periods  of  time. 

It  was  then  requested  that  the  necessary  modifications  be  included  to  take 
into  account  the  drift  of  mean-geostationary  satellites  with  time  for  periods 
of  up  to  several  months.  Positional  accuracy  of  sub-satellite  longitude  of 
approximately  one  degree  was  deemed  sufficient.  However,  satellite  latitude 
positional  accuracy  of  .2  degrees  was  necessary  because  the  conversion  tables 
are  generally  more  dependent  upon  satellite  latitudes. 

Rather  than  make  the  corrections  to  the  original  LOKANGL  program, 
it  was  decided  that  it  would  be  easier  and  more  efficient  to  modify  a version 
of  LOKANGL  called  ROPLOK . (Rapid  Orbit  Prediction  Table).  (For  a further 
description  of  ROPLOK  refer  to  references  1,  2,  3,  at  the  end  of  this  section.) 

ANALYSIS  AND  RESULTS 

After  the  necessary  changes  were  made  to  ROPLOK,  it  was  combined  with 
LOKANG2.  The  ROPIOK  program  also  contained  unnecessary  and  duplicate 
statements  which  were  eliminated.  The  input  to  this  program  was  changed  to 
accept  a specific  5-card  input  deck.  It  was  noted  that  both  programs  originally 
had  some  erroneous  statements  which  were  corrected  and  tested  many  times  to 
ensure  the  accuracy  of  the  algorithm.  To  do  this,  the  program  was  run  on 
different  input  samples  and  cross-checked  with  theoretical  results.  The  new 
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combined  program  produces  a Rapid  Orbit:  Prediction  Table,  a table  of  M-values, 
look  angles  for  a geostationary  satellite  as  a function  of  satellite  latitude 
and  longitude,  and  a Faraday-to-TEC  conversion  table. 

All  subroutines  in  ROPLOK  were  modified  to  minimize  compile  and 
execution  time  as  well  as  storage  space.  The  unnecessary  statements  and 
unused  variable  were  deleted.  Nine  subroutines  from  L0KANG2  were  added  to 
ROPLOK  program.  These  subroutines  compute  the  M-value  table.  The  resulting 
program  is  called  R0PLK2. 

In  order  to  run  program  R0PLK2  on  UNIVAC  1108  machine,  similar 
changes  made  to  LOKANG2  were  made  to  R0PLK2,  and  the  program  was  con- 
verted from  029  keypunch  code  to  026.  Again  all  I/O  is  done  via  GREAD 
and  GPRINT  subroutines.  These  are  the  UNIVAC  1108  I/O  routines  at 
AWS-AFGWC. 

The  following  M-Value  routines  from  LOKANG2  were  added  to  ROPLOK 

program: 

1.  MGEOS  ( MCON ) 

2.  GEO 

3.  EIELDG 

4.  FIELD 

5.  HOUR 

6.  HEED 

7 . TABLE 

8 . XYZ 

9.  AZIM 

In  order  to  optimize  the  program,  the  following  ROPLOK  routines 
were  deleted  or  altered  as  follows: 
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A 


SUBROUTINE 

ACTION 

REASON 

PRTCON 

Deleted 

This  ubroutine  only  prints  the  J(L,M) 
Coeff.,  LAMDA  ( L , M)  Coeff.,  Subscripts 
(6  Sets),  other  constants,  and 
Perturbations  Present. 

LABEL 

Deleted 

Called  only  by  PRTCON  (Deleted  Above) 

INPREF 

Deleted 

Used  for  Tape  7 (Name  List)  not  desired 

MEAN IN 

Altered 

Several  lines  were  deleted  since 
they  checked  for  U or  X (more  than 
one  element  set) . 

POVLIN 

Deleted 

Called  only  by  MEANIN  if  more  than 
one  element  set. 

CCDAM 

Deleted 

Reads  in  a second  set  of  elements. 

PVTELM 

Deleted 

Called  by  POVLIN  and  CCDAM  (Both 
deleted) . 

NODSTR 

Deleted 

Called  only  by  MOON  and  is  therefore 
incorporated  therein. 

XCROSS 

Deleted 

Called  only  by  EPHSET  and  is  therefore 
incorporated  therein. 

BLOCK  DATA 

Deleted 

Inserted  into  main  program. 

The  final  step  was  to  guarantee  that  the  program  be  compatible 
with  the  UNIVAC  1108  at  the  AWS-AFGWC.  It  was  determined  that  an  1108 
compiler  or  simulated  compiler  did  not  exist  at  AFGL.  Several  routines 
which  simulate  1108  input/output  functions  did  exist.  However,  the 
primary  one  of  interest,  XOWAIT,  was  not  running  properly  on  the  AFGL  6600 
system.  It  was  then  decided  to  perform  all  input/output  operations  by 
using  BUFFER  IN,  BUFFER  OUT,  ENCODE  and  DECODE  routines  (i.e.,  GREAD  and 
GPRINT  subroutines) . Special  routines  were  designed  to  simulate  the 
UNIVAC  1108  input/output  system  routines.  One  last  problem  had  remained, 
that  of  temporary  mass  storage  units.  This  problem  was  overcome  by  partial 
storage  of  the  values  in  the  central  memory  and  re-allocating  the  storage 
as  it  was  needed. 
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14. 


Persistence,  Recurrence  And  Runs 
Of  Various  Weather  Elements 


INITIATOR:  I.  LUND 

PROBLEM  NO:  4841  PROJECT  NO:  8624 


BACKGROUND 

Past  climatological  studies  have  been  performed  to  investigate 
various  weather  elements.  A software  system  has  been  available  to  extract 
the  dat  and  perform  the  calculations.  Portions  of  the  software  have  been 
determined  to  be  both  slow  and  cumbersome  to  use.  New,  more  efficient 
algorithms  have  been  added  and  program  code  improvements  have  been  made 
to  part  of  the  software  system.  Additional  effort  is  required  to  complete 
•■he  system  optimization. 

ANALYSIS  AND  RESULTS 

The  data  base  for  past  climatological  studies  consisted  of  13 
years  of  hourly  weather  observations  from  nine  weather  stations  on  the 
east  coast  of  the  United  States.  The  stations  are: 

1.  Andrews  AFB , D.C. 

2.  Raleigh-Durham,  NC  (Airport) 

3.  Philadelphia,  PA  (International  Airport) 

4.  Byrd  Field,  Richmond,  VA 

5.  Washington,  D.C.  (National  Airport) 

6.  New  York  (LaGuardia  Airport) 


■.»  •ii-v-m*  r-n 
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7.  Newark,  NJ  (Airport) 


B.  Baltimore,  MD  (Friendship  International  Airport) 

9.  New  York  (J.F.K.  International  Airport) 

Included  in  the  computations  for  these  studies  was  the  computation  of 
frequency  and  cumulative  frequency  distributions  for  various  weather 
elements  for  each  station  iind  all  combinations  of  stations.  The  prob- 
abilities of  occurrence,  persistence  and  recurrence  of  specified  ranges 
of  values  for  different  seasons  were  calculated.  Based  on  these,  prob- 
abilities and  distributions  of  various  weather  elements  can  be  categorized. 

Weather  elements  of  interest  include; 

• Ceiling 

• Visibility 

• Wind  Speed 

• Temperature  (Dry  Bulb) 

• Dew  Point 

• Sea  Level  Pressure 

• Total  Cloud  Amount 

• Opaque  Cloud  Amount 

Previous  analyses  produced  a set  of  programs  which  were  used  to 
organize  the  data  and  to  calculate  the  probabilities.  From  examination 
of  the  software,  it  was  easy  to  identify  areas  which  could  be  improved  to 
facilitate  the  analysis  of  the  climatological  data.  Two  of  the  programs 
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in  thi>  analysis  software  package  have  already  been  updated.  it  is  antici- 
pated that  future  similar  work  will  improve  the  total  software  and,  with 
tlu>  more  efficient  procedures,  facilitate  the  analysis.  Specifically,  the 
software  has  been  improved  in  four  areas: 

• Computer  Run  Time 

• Flexibility 

• Code 

• Output 

The  original  algorithm  used  to  calculate  the  joint  probability 
distribution  was  slow.  For  a sample  case,  the  execution  time  was  950 
seconds.  The  new  algorithm  performs  the  calculations  in  56  seconds. 


In  a binary  system  (that  is,  either  the  event  occurred  or  it  didn't) 

we  have 


p(X  ) 


N 

Total  number  of  points  in  data  set 


p(X  = X. ) 

1 


where  N = number  of  occurrences  for  an  event  (number  of  favorable  cases) 
X = random  variable  in  the  set  of  observation 


To  simplify,  assume  X = 1 for  an  event  occurring,  then  p(X. ) = p(l  = X^) . 

For  a two  dimensional  case,  we  have 

pfX^,  X^)  = p(X  = and  X = X^)  = pfXjApfX.) 

= p (X , = 1 and  X . = 1 ) 
l 1 

For  a three  dimensional  case«  we  have 

p(X.,  X.,  X)  = p (X , = 1 and  X,  = 1 and  X = 1) 

1 ) K 1 ] K 

= p (X . ) Ap  (x  . ) Ap  (X  ) 
i j k 
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number  of  times  X.,  X.  and  X.  all  occur  at  the  same  time 

i J k 

Total  number  of  observations 


Similar  equations  can  be  written  for  four  through  nine  dimensional  cases. 


For  this  study,  there  was  data  from  nine  locations.  Probabilities 
were  desired  for  each  site  and  all  combinations  of  sites.  That  is,  the 
probability  of  occurrence  at  sites  up  to  nine  at  a time  were  required.  Two 
basic  counters  were  required:  one  to  indicate  number  of  data  points 

available  for  each  site  and  each  combination  of  sites;  the  other  was  a 
counter  of  occurrence,  i.e.,  number  of  times  site  and/or  sites  met  the 
criterion,  for  each  combination  of  sites.  The  total  number  of  counters 
required  is 


2 


= 1022 


Originally  the  program  was  written  to  search  through  all  combina- 
tion of  sites  to  determine  if  the  criteria  was  met  for  that  combination. 
However,  if  for  example  only  two  sites  met  the  criterion  there  is  no  reason 
to  search  through  the  combinations  of  three  or  more  sites.  This  portion, 
i.e.,  the  section  to  determine  how  many  and  which  sites  had  an  occurrence  , 
was  rewritten.  If  no  sites  indicated  an  occurrence,  no  searching  would  be 
performed.  The  total  possible  data  counters  would  be  increased  and  addi- 
tional data  would  be  examined.  Second  if  all  sites  had  an  occurrence,  all 
counters  would  be  increased.  Third,  combinations  would  be  searched  up  to 
the  number  of  combinations  possible.  That  is,  if  for  example  three  sites 
indicated  an  occurrence,  no  combination  of  four  or  more  would  be  searched. 
Further,  sites  were  ordered  so  that  additional  searching  would  be  minimized. 
For  example,  assume  sites  2,  3,  4 and  7 indicated  an  occurrence.  The  counter 
array  would  be  increased  based  on  an  algorithm  keying  on  the  site  iden- 
tifying number.  Further  for  the  cases  where  the  number  of  sites  at  a 
particular  time  was  five  through  eight,  keying  on  the  sites  without  an 
occurrence  could  be  performed.  That  is,  increase  all  counters  by  one, 
then  decrease  the  counters  for  the  sites  where  there  was  no  occurrence. 

Hence,  if  eight  sites  reported  an  occurrence,  increase  all  counters  then 
decrease  all  counters  including  M,  the  site  with  no  occurrence.  Inclusion  of 
these  simple  changes  produced  a significant  decrease  in  execution  time. 
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Not  all  changes  are  expected  to  produce  such  a dramatic  effect; 
however,  additional  savings  are  anticipated.  The  original  programs  required 
the  repunching  and  substitution  of  many  cards  within  the  program  to  change 
the  running  of  a category  (range,  season,  element  combination)  of  interest. 
The  programs  were  changed  to  allow  for  various  categories  to  be  run  with 
minor  changes  in  the  input.  The  program  code  was  cumbersome  and  at  times 
redundant.  Portions  of  the  code  were  rewritten,  thereby  reducing  code  and 
improving  the  logic  of  the  original  program.  Various  additional  output 
parameters  were  added  to  the  program.  These  parameters,  such  as  number  of 
occurrences,  provided  more  input  to  the  categorization  procedure. 

No  major  difficulties  are  anticipated  in  the  continuation  of  this 
effort.  The  continued  effort  will  include  the  following  steps:  review  the 

total  system,  data  base  and  software;  in  conjunction  with  analyst  input, 
determine  the  amount  of  use  this  system  will  be  exercised  for  future 
analyses;  then,  depending  upon  anticipated  activity,  either 

1)  optimize  the  data  base  handling, 

or  2)  continue  to  optimize  individual  program, 

or  3)  combine  portions  of  individual  programs  to 
optimize  the  system, 

or  4)  combinations  of  1),  2)  and  3). 

That  is,  some  efficient  procedures  for  calculating  probabilities  of 
occurrence,  persistence  and  recurrence  of  specified  ranges  have  already 
been  identified  and  included.  Additional  improvements  could  be 
identified  and  added  to  the  total  system  to  enhance  its  performance  in 
providing  a detailed  analysis  of  the  climatological  data. 
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15. 


White  Sands  Missile  3 x 80  Range  Tape  Conversion 


INITIATOR:  E.  ROBINSON 


PROBLEM  NO:  4904  PROJECT  NO:  0001 


BACKGROUND 

Data  tapes  were  received  from  White  Sands  Missile  Range  (WSMR) . 
These  tapes,  generated  on  a tJNIVAC  1108,  contain  multiple  radar  station 
data.  Each  record  is  composed  of  80  words  and  each  block  of  data  con- 
tains 3 records.  Thus,  if  there  are  2 radar  stations.  Record  one  will 
contain  the  header  for  Station  #1;  Record  two,  the  header  for  Station  “2; 
Record  three,  data  for  Station  #1:  Record  four,  data  for  Station  #2;  etc. 

ANALYSIS  AND  RESULTS 

A program  was  supplied  to  read  tapes  containing  data  for  a single 
radar  station.  This  program  was  modified  to  accept  the  additional  radar 
stations.  The  program  originally  buffered  in  1 block  of  input  data  and 
set  up  a loop  to  convert  the  36  bit  floating  point  words  into  60  bit 
floating  point  words  in  a manner  similar  to  that  explained  in  Section  5 
of  this  report.  Once  a block  of  data  was  processed,  the  data  was  output 
to  a tape  file  and  to  the  printer;  whereupon  the  program  returns  to  read 
more  data. 

The  program  was  modified  by  first  reading  from  cards  the  number 
of  radar  stations  contained  on  the  tape,  and  then  setting  up  a 3-dimensional 
array  for  output  instead  of  the  original  6 x 16  array.  This  array  is  in 
the  form  of  (I  x 6 x 16) . I corresponds  to  the  number  of  Radar  Stations. 
Each  input  block  contains  3 data  records.  The  data  is  made  available  for 
output  to  separate  files  for  each  Radar  Station  only  after  the  array  is 
filled  with  data.  The  array  is  then  cleared  after  output  and  new  data 
processed . 
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The  coordinates  of  the  origin  and  of  the  radars  were  extracted 
from  the  headers.  This  information  was  not  contained  in  the  original 
program.  A conversion  for  integer  numbers  was  included  in  the  conversion 
package  to  obtain  these  data  values.  This  conversion  method  is  also 
explained  in  Section  5. 

Five  tapes  were  processed.  These  tapes  included  both  single 
and  multiple  radar  stations.  Problems  with  the  units  of  the  time 
value  associated  with  each  data  point  were  encountered.  A time  conversion 
algorithm  was  included  to  obtain  the  proper  units. 
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16. 


Coordinate  Conversions  And  Graplical  Presentations 


INITIATOR:  W.  BARRON 

PROBLEM  NO:  4908  PROJECT  NO:  433L 


BACKGROUND 

Often  computer-based  graphical  presentations  are  useful  for  analysis 
and  as  pictorial  descriptions  of  systems  under  study.  They  are  used  to  dis- 
play values  commonly  associated  with  different  but  equivalent  geometric 
bases.  For  example,  there  exists  a requirement  to  display  the  astronomical 
altitude-azimuth  coordinate  system  and  the  astronomical  hour  angle-declina- 
tion coordinate  system,  centered  around  any  specific  geographic  location. 

On  the  display  or  graph,  an  example  of  which  is  shown  in  Figure  16-1,  various 
other  information  such  as  satellite  tracks  can  be  displayed  to  aid  in  the 
analyses.  The  algorithm  performing  the  coordinate  conversion  is  useful  to 
locate  accurately  any  object  in  either  coordinate  system  given  one  set  of 
coordinates . 

The  celestial  sphere  provides  a convenient  framework  for  fixing  r':.- 
relative  positions  of  the  heavenly  bodies.  A point  on  the  celestial  sphere 
directly  above  the  observer's  position  is  the  zenith.  The  point  on  the 
celestial  sphere  directly  opposite  the  zenith  is  the  nadir.  The  plane 
perpendicular  to  zenith-nadir  line  at  the  observer's  position  intersects 
the  celestial  sphere  in  the  horizon.  This  is  the  great  circle  whose  points 
are  90  degrees  from  the  zenith  everywhere.  A northward  prolongation  of  the 
earth's  axis  of  rotation  pierces  the  celestial  sphere  at  a point  called  the 
north  celestial  pole.  Tt  is  marked  by  a moderately  bright  star  --  Polaris. 
However,  there  is  no  pole  star  in  the  southern  hemisphere.  The  observer's 
meridian  is  a great  circle  through  the  zenith  and  the  north/south  poles. 

The  position  of  a celestial  object  is  given  by  its  altitude  and  its  azimuth. 
The  latitude  is  the  angle  in  degrees  measured  from  the  horizon  along  the 
vertical  circle  to  the  celestial  object.  The  zenith  distance  is  90  degrees 
minus  latitude.  The  azimuth  is  the  angle  in  degrees  measured  from  the 
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north  point  toward  the  east  along  the  horizon  to  the  foot  of  the  vertical 
circle  passing  through  the  celestial  object.  The  great  circle  which  passes 
through  the  north  pole  and  a celestial  object  and  is  perpendicular  to  the 
equator  is  called  an  hour  circle.  The  declination  of  this  celestial  object 
is  the  angular  distance  in  degrees  from  the  equator  along  the  hour  circle 
to  the  celestial  object. 

In  order  to  be  able  to  describe  the  position  of  a heavenly  body 
in  the  sky,  it  is  convenient  to  suppose  the  inner  surface  of  the  celestial 
sphere  to  be  marked  off  by  imaginary  circles  (e.g.,  the  equator),  the 
meridians  and  the  parallels  of  latitude  which  serve  as  a system  of  reference 
for  the  latitude  and  longitude  of  a point  on  the  surface  of  the  earth. 

Several  distinct  systems  of  such  circles  are  used  in  astronomy. 

The  conversion  from  one  coordinate  system  to  another  is  done  by  some  simple 
equations.  The  primary  coordinate  conversion  of  interest  is  the  altitude- 

azimuth  coordinate  system  to  the  hour  angle-declination  coordinate  system 
and  vice  versa. 

The  objective  of  this  problem  addressed  under  Contract  No.  F19628- 
76-C-0203  is  two  fold:  first,  to  write  a program  for  the  conversion  of  the 

Astronomical  Altitude-Azimuth  Coordinate  System  to  the  Astronomical  Hour 
Angle-Declination  Coordinate  System,  and  vice  versa,  for  any  geographic 
location.  The  output  of  this  computation  is  to  be  used  i,  aiding  the  in- 
stallation and  operation  of  radio  telescope  systems  to  be  located  through- 
out the  world.  Second,  write  plotting  software  to  display  the  two  coordinate 
systems  centered  about  a given  geographic  latitude. 
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ANALYSIS  AND  RESULTS 


Based  on  the  geometry  shown  in  Figure  16-2,  the  equations  used 
in  the  coordinate  conversions  are  derived  as  follows: 


Figure  16-2  Geometry  for  Coordinate  Conversion 


Equation  Is  azimuth-zenith  angle  to  hour  angle-declination 


For  a given 


c = co-latitude 
A = azimuth 
b = zenith  angle. 


use  the  following  equations 


cos  a = cos  b cos  c + sin  b sin  c cos  A 


„ cos  b - cos  a cos  c 

cos  B = ; 

sin  a sin  c 


to  compute 

B =,  hour  angle 
6=9)-  a = declination 
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Equation  2:  hour  angle-declination  to  azimuth- zenith  angle 


For  a given 

c = co-latitude 

B = hour  angle 

a = 90  - 6 = co-declination, 


use  the  following  equations 


cos  b = cos  a cos  c + sin  a sin  c cos  B 


cos  A 


cos  a - cos  b cos  c 
sin  b sin  c 


to  compute 


A = azimuth 
b = zenith  angle 


The  coordinate  conversion  program  is  exercised  to  produce  the  plot 
of  the  Astronomical  Altitude-Azimuth  Coordination  System  and  the  Astronomical 
Hour  Angle-Declination  Coordinate  System  centered  about  any  given  latitude 
including  zero.  This  plot  may  be  displayed  on  30- inch  paper  as  well  as 
12-inch  paper  with  multicolor  options. 


The  software  is  written  so  that  modifications  and  extensions  can 
readily  be  made  for  future  use.  Preliminary  discussions  have  focused  on 
the  following  expansions  of  the  capability  of  the  software: 


1)  Display  a satellite  track  (centered  about  some  geo- 
graphic latitude)  both  on  the  plot  and  in  digital 
form. 

2)  Magnify  portions  of  any  plot  (for  a given  altitude, 
declination,  and  hour  angle). 

3)  For  a given  time  range  (usually  a month) 

and  latitude,  plot  the  visible  constellation. 

4)  Add  additional  coordinate  system  to  the 
software  for  conversion  and  display  information; 
e.g..  Right  Ascension  - Declination 

Latitude  - Longitude 

Galactic  Latitude  - Galactic  Longitude 
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17. 


Contour  Plotting  Routines 


INITIATOR:  E.  SHETTLE 

PROBLEM  NO:  4909  PROJECT  NO:  7621 


BACKGROUND 

A package  of  subroutines  was  developed  to  plot  a contour  map 
for  a given  function.  The  main  subroutine,  CONTOUR,  accepts  as  input 
an  array  of  function  values  f(X,  Y)  on  a grid  of  points  (X(I),  Y(J)) 
and  plots  a contour  map  (or  maps)  of  the  function.  Two  interpolation 
algorithms  are  provided  to  find  values  of  f(X,Y)  when  X and  Y do  not 
fall  on  the  specified  data  points.  More  interpolation  procedures  may  easily 
be  added  to  the  existing  ones.  For  eight  X values  and  nine  Y values  (i.e., 
72  f(X,  Y)  values),  this  subroutine  produces  two  sets  of  contour  plots,  one 
with  two  contour  values  and  the  other  with  seven  contour  values,  in  five 
CPU  seconds. 

ANALYSIS  AND  RESULTS 

This  subroutine  package  is  composed  of  a drive  subroutine  and 
nine  other  subroutines.  The  driver  subroutine  accepts  the  array  values 
as  well  as  a set  of  parameters  which  determines  type  of  interpolation 
desired,  and  the  scale  and  level  of  contours  to  be  displayed.  The  output 
is  the  plot  of  a contour  map  of  a given  function.  The  function  of  each 
subroutine  is  to  evaluate  a set  of  variables  and/or  to  perform  a chain  of 
interactive  operations. 

Individual  subroutines  were  written  to 

1 . Accept  a set  of  parameters  and  check  for 
possible  errors 

2.  braw  frame,  axes  and  labels 
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3.  Label  contours 

4.  Generate  a set  of  contour  values 

5.  Generate  second  set  of  contour  values 

6.  Interpolate  the  values  of  the  function  when  X 
and  Y do  not  fall  on  the  specific  data  points 
for  either  linear  or  logarithmic  values 

7.  Check  continuity  of  the  contour  curves 

8.  Relocate  coordinates  of  contour  map  to  smooth 
the  curve 

9.  Sort  values  and  obtain  least  upper  bound  and 
greatest  lower  bound  for  a specified  value 

The  interpolation  algorithms  used  to  find  values  of  f(X,  Y) , 
when  X and  Y do  not  fall  on  the  specified  data  points  X^,  (i  = 1,  . ..,  NX) 

and  Yj , (j  = 1,  2,  NY),  are  given  below. 

When  f is  linear  in  both  X and  Y the  following  equations  are  used: 


f(X,  Y)  = (1  - h ) (1  - h ) f.  . + (1  - h ) h f.  . 

x y l,  D x y i,  :+l 


+ h (1  - h ) f.  . + h h f. 

x y 1+1,3  xyi+1,3+1 


where  h = (x-x.)/(x.  , -x.) 

x l l + 1 1 


hy  = (y  - YjJ/tYj  + ! - Vj) 

When  log  (f)  is  linear  in  both  X and  Y,  the  following  equations  are  used: 

Log  f(x,  y)  = (1  - h ) (1-h)  Log  f.  . + (1  - h )h  Log  f. 

1 x y 1,  : x y ^ i,  3 + 

+ h (1  - h ) Log  f.  , . + h h Log  f.  . . . , 

x y ^i+l,  3 xy  l +1,  3+1 


where  h and  h are  the  same  as  above, 
x y 
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This  subroutine  package  was  designed  to  be  user  oriented.  Much 
care  was  taken  in  the  writing  of  the  software  to  facilitate  its 
use.  The  software  code  itself  has  much  documentation  via  comment  cards. 
The  following  was  taken  from  the  program  documentation.  It  is  a detailed 
explanation  of  the  input  required.  It  also  indicates  that  this  package 
is  sufficiently  general  in  nature  to  have  varied  applications  and  that 
with  minor  modifications  the  range  of  applications  could  be  expanded 
further . 


The  input  to  this  subroutine  is  provided  by  the  user.  The  principal 
subroutine  CONTOUR  is  independent  of  the  main  program  and  has  no  READ  or  WRITE 
statements.  It  is  the  user's  responsibility  to  provide  the  input  to  this  sub- 
routine by  passing  it  as  the  parameters  when  a call  is  made  to  this  subroutine. 
The  input  to  this  subroutine  is  explained  in  details  on  the  following  pages. 

For  simplicity,  the  following  codes  are  used  for  variable  des- 
cription : 


RU  = Returned  Unchanged 
R = Real  Value 
INT  = Integer  Value 

Call  to  this  subroutine  is  made  by: 


CALL  CONTOUR  (F,  X,  Y,  NX,  NY,  INTERP,  NFC,  FCON,  XCALE, 
XNAME,  YNAME , HEADER,  INFO,  XCONTR,  YCONTR,  BEGINX, 
BEGINY,  MESS) 


Variable  Description: 

F = A2  - Dimensional  array  such  that: 

F (I,  J)  = F(Xi(  Y.)  with  X.  = X ( I ) and  ¥ ^ = Y(J) 

Dimensions  of  this  matrix  is  NX  by  NY. 

For  every  call,  the  dimension  of  F in  the  calling 
program  must  be  exactly  F(NX,  NY).  (RU)  (R) 
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X = A vector  of  X values  sorted  in  ascending  or 

descending  order  for  which  F(X,  Y)  is  specified. 
Dimension  of  X is  NX,  (RU)  (R) 

Y = The  corresponding  vector  of  Y values  where 
F(X,  Y)  is  given.  Dimension  of  Y is  NY. 

(RU)  (R) 

NX  = Number  of  points  in  the  X vector  (RU)  (INT) 

NY  = Number  of  points  in  the  Y vector  (RU)  (INT) 

INTERP  = The  name  of  an  external  subroutine  to  be 

used  for  interpolation  when  the  coordinates 
of  a contour  curve  do  not  fall  on  the  spe- 
cified data  points.  (RU) 


= INTRP0,  if  f is  linear  in  both  X and  Y. 


f(X,  Y)  = (1  = h ) (1  - h ) f.  . + (1  - h ) h f. 


y i-  i 


x y i,  j + 1 


+h  (1  - h ) f.  . ,+hh  f . , . 

x y i+I,  3 xy  l+l,  3+1 


whe re:  h = ( x - x . ) / ( x . , . -x.) 

X 11  + 1 1 


h = (y  - y . )/(y . , , - y . ) 

y ']  1+1  1 


= INTRP1 , if  bog  (f)  is  linear  in  both  x and  y. 


Log  f (x,  y)  = (1  - h ) (1  - h ) Log  f. 

x y i,j 


+ (1  - h ) h Logf. 

x y l,  ] + 1 

+ h (1  - h ) Log  f. 

x y 3 l + 1,  j 


+ h h Log  f . , . . 

xy  l+l, 1+1 


where:  h and  h are  the  same  as  above, 

x y 

These  subroutines  are  provided  for  the  use  in  this  program.  The  user  must 
pass  one  of  the  subroutine  names  to  this  program  and  it  must  be  declared 
external  in  the  calling  program. 
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NFC  = Number  of  contours  to  be  plotted.  (INT) 


= Positive,  NFC  contours  at  the  values  given  by  FCON(I),  1=1,  2, 

NFC  are  drawn.  (RU) 

= 0,  the  contours  such  that  FC0N(1)  + N*FC0N(2)  = FC0N(3)  where 

N = 0,  1,  2,  ....  are  drawn.  Upon  return  NFC  = Number  of  contours. 

= Negative,  -NFC  contours  between  FMIN  and  FMAX,  inclusive,  where  FMIN 
and  FMAX  are  the  smallest  and  the  largest  values  in  the  array  F(I,  J) 
are  drawn.  Upon  return  NFC  = ABS(NFC). 

FCON  = An  array  of  contour  values.  (R) 

a - if  NFC  > 0,  then  FCON  must  contain  NFC  contour  values.  (RU) 

b - if  NFC  = 0,  then  FC0N(1)  = Initial  contour  value,  FCON (2)  = Increment 
value  and  FCON(3)  = Maximum  contour  value.  Upon  return  FCON  contains 
the  evaluated  contour  values. 

c - if  NFC  < 0,  then  FCON  is  evaluated  between  FMIN  and  FMAX.  If  NFC  = -1, 
then  FCON ( 1 ) = FMIN.  Upon  return  FCON  contains  the  evaluated  contour 
values . 

Note:  Dimension  of  FCON  is  L,  where  L is  as  follows: 

L = NFC  if  NFC  > 0. 

L = !NFC  j if  NFC  < 0. 

L = INT ( FCON (3)  - FCON ( 1 ) ) /FCON ( 2)  + 1 if  NFC  = 0. 

XSCALE  (I)  = A vector  defining  the  X-axis  of  the  plot.  (RU) (R) 

(1)  = Left  most  X value  to  be  written  on  the  X-axis  of  the  plot. 

(2)  = Right  most  X value  to  be  written  on  the  X-axis  of  the  plot. 

(3)  = Number  of  divisions  per  10  inches 
= 10.  means  1 tick  mark/inch. 

= 5.  means  1 tick  mark/2inches . 

= 20.  means  1 tick  mark/half  inch. 

(4)  = Length  of  X-axis  in  inches 

Note:  Dimension  of  XSCALE  is  4. 
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YSCALE  (I)  = A vector  defining  the  Y-axis  of  the  plot.  (RU) (R) 

(1)  = Bottom  most  Y value  to  be  written  on  the  Y-axis  of  the  plot. 

(2)  = Top  most  Y value  to  be  written  on  the  Y-axis  of  the  plot. 

(3)  = Same  as  XSCALE(3). 

(4)  = Length  of  Y-axis  in  inches. 

£ 10  for  12-inch  width  paper. 

< 28  for  30-inch  width  paper. 


XNAME  = An  alphanumeric  array  containing  up  to  40  characters  to  label 
the  X-axis. 

Note:  Dimension  of  XNAME  is  4. 

YNAME  = An  alphanumeric  array  containing  up  to  40  characters  to  label 
the  Y-axis. 

Note:  Dimension  of  YNAME  is  4. 

HEADER  (I)  = An  alphanumeric  array  containing  up  to  110  characters.  First 
3 words  are  used  as  identification  for  plot.  Words  4 through 
11  are  printed  above  the  contour  plot. 

(1)  = 10H  ...  Programmer's  name. 

(2)  = 10H  ...  Problem  number. 

(3)  = 10H  ...  Information  meaningful  to  the  programmer. 

Note:  Not  all  words  from  (4)  to  (11)  must  necessarily  be  used. 

Dimension  of  HEADER  is  11. 

INFO  (I)  = Control  Parameters.  (RU)  (int) 

(1)  = Number  of  characters  in  XNAME. 

(2)  = Number  of  characters  in  YNAME. 

(3)  = Number  of  characters  in  HEADER (4 ) through  HEADER (N) , where 

4 < N < 11. 

(4)  = Control  for  using  symbols. 

= 0,  line  plot  without  symbols. 

= 1,  line  plot  with  a symbol  at  every  point. 

= Negative  will  suppress  lines  between  the  points. 


R2 


(5)  = Symbol  code  used  for  plotting. 


= 3 — symbol  = +. 

=4  — symbol  = X. 

= 11  — symbol  = * 

(6)  = Number  of  times  this  subroutine  will  be  called.  This  determines 

the  number  of  frames. 

(7)  = Total  length  of  X-axis  in  inches  for  all  frames.  That  is,  sum 

of  all  XSCALE(4)  for  all  calls. 

Note:  INF0(6)  and  INF0(7)  must  be  the  same  for  all  calls. 

Dimension  of  INFO  is  7. 

XCONTR  = Dummy  array  to  store  X values  of  the  same  level  contour  curves.  (R) 

YCONTR  = Dummy  array  to  store  the  corresponding  Y values  of  the  contour 
curves.  (R) 

BEGINX  = Dummy  array  to  store  the  starting  X coordinates  of  contour  curves 
of  the  same  level.  (R) 

BEGINY  = Dummy  array  to  store  starting  Y coordinates  which  correspond  to 
X coordinates  in  BEGINX.  (R) 


Note:  Dimensions  of  these  4 dummy  arrays  are  to  be  determined 

by  the  user. 


MESS  (1)  = Error  or 

warning  check.. 

(INT) 

(2)  = Dimension 

of 

XCONTR 

and 

YCONTR 

declared 

in  the 

calling 

program 

(3)  = Dimension 

of 

BEGINX 

and 

BEGINY 

declared 

in  the 

cal 1 ing 

program 

IMPORTANT  POINTS 

1.  Plotter  requires  the  information  in  HEADER ( 1 ) through 
HEADER(3).  It  must  be  the  same  for  all  calls. 

2.  This  subroutine  may  be  called  as  many  times  as  desired.  Define 
the  parameters  as  described  above  for  every  call.  For  every 
call  the  dimension  of  F in  the  calling  program  must  be  exactly 
NX  by  NY. 

3.  CALL  ENDPLT  must  be  used  to  obtain  a plot.  Insert  this  state- 
ment before  the  end  of  the  main  program. 

4.  F values  must  be  passed  as  follows: 
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X (1 ) 

X (2 ) 

— 

X (NX) 

Y(l) 

F(l, 

1) 

F(2, 

1) 

— 

F (NX , 

1) 

Y (2 ) 

F (1 , 

2) 

F (2 , 

2) 

— 

F (NX , 

2) 

Y (NY) 

F (1 , 

NY) 

F (2  , 

NY) 

— 

F (NX , 

NY) 

The  following  read  statement  would  do  the  job: 

READ  100,  ( ( F ( I , J) , 1=1,  NX),  J = 1,  NY) 

100  FORMAT  (nx  Fw.d) 
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18. 


Program  And  Tape  Conversion  For  ATS-G  Satellite  Data 


INITIATOR:  D.  HARDY 

PROBLEM  NO:  4923  PROJECT  NO:  7661 


BACKGROUND 

Tapes  of  the  data  generated  by  the  ion  and  electron  spectrometers 
on  board  the  ATS-6  satellite  were  furnished.  These  tapes  were  generated  on 
the  UNIVAC  1108  computer  which  has  a 36-bit  word  size.  A program  that  reads 
the  data  tapes  and  obtains  particle  counts,  velocities,  and  distribution 
functions  was  obtained  from  University  of  California,  San  Diego  (UCSD) . 

This  program  was  prepared  on  a CDC  3600  computer  with  a 48-bit  word.  It 
was  desirable  to  have  the  program  operational  on  the  CDC  6600  computer, 
and  to  be  able  to  read  the  tapes  with  this  program. 

ANALYSIS  AND  RESULTS 

This  task  is  two-fold  in  nature,  since  neither  the  tape  nor  the 
program  was  compatible  with  the  CDC  6600.  Modifications  were  made  to 
the  program,  and  then  subroutines  were  added  to  the  program  to  read  the  tape 
properly . 


The  program  was  first  duplicated  and  compiled.  All  compilation 
errors  were  removed.  Logical  operations  had  to  be  checked  and  modified 
as  needed.  The  output  header  formatting  had  to  be  changed  because  of 
differences  in  the  word  sizes  for  the  two  computers.  A routine  for  tape 
positioning  was  written  to  replace  an  internal  reference  to  tape  position- 
ing available  on  the  CDC  3600,  but  not  available  on  the  CDC  6600. 

The  major  problems  encountered  arose  when  processing  the  first 
records  of  the  tape  began.  Each  file  is  preceded  by  a record  in  external 
Binary  Coded  Decimal  (BCD) . Parity  cannot  be  changed  during  the  execution 
of  a program.  Therefore,  for  the  processing  of  input,  a conversion  had  to 
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be  written  to  change  the  characters  from  external  BCD  to  internal  BCD,  with 
minor  modifications  for  masking  and  decoding  of  the  record.  The  starting 
time  is  obtained  from  this  record;  thus,  it  was  an  essential  step  in  the 
processing. 

In  addition  to  this  file  header  record  problem,  the  data  records 
caused  problems  because  of  lack  of  compatibility  of  word  sizes.  The  pro- 
gram supplied  assumed  a tape  word  size  of  36  bits,  but  assumed  an  internal 
computer  word  size  of  48  bits.  However,  60  bits  are  normally  used  when 
programming  on  a CDC  6600.  Therefore,  we  chose  to  unpack  48  bits  of  tape 
information  into  a 60-bit  word,  right  justified,  to  obtain  a compatible 
process,  and  let  the  program  conversion  do  the  rest.  This  was  considered 
preferable  to  modifying  each  step  of  the  program  to  account  for  the  extra 
12  bits  in  the  CDC  6600  word  size.  This  operation  was  accomplished  by 
setting  up  loops  and  shifting  three  60  bit  words  of  information  into  five 
60  bit  words  containing  48  bits  of  information,  right  justified,  with  zero 
fill  for  the  remainder  of  the  word. 

The  initial  output  values  obtained  were  relatively  smaller  than 
expected,  and  further  inspection  of  the  program  showed  that  there  is  an 
option  to  obtain  the  counting  rates  either  in  counts/second  or  Log^  (counts/ 
sec) . The  sample  listing  supplied  was  in  counts/second.  The  Log  conversion 
was  by-passed,  and  a duplicate  of  the  sample  listing  was  obtained.  Other 
periods  of  interest  were  processed  with  satisfaction  and  duplicate  copies  of 
the  new  deck  were  made  for  future  use. 

In  general,  this  program  reads  data  generated  from  the  University 
of  California,  San  Diego  (UCSD)  experiments.  Observations  are  made  of 
trapped  particles  and  their  interrelations  in  the  following  regions: 

• trapped  radiation  belts  toward  the  equator. 

• N-S  auroral  zones 

• Daytime  magnetopause  and  interplanetary  sphere 

• Nighttime  tail  plasma  sheets. 
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Also  observed  in  these  regions  will  be  the  energy  spectrum  of  ions  above 
50  eV,  density  versus  magnetic  field  depression,  and  fluxes  of  auroral 
particles  and  their  associations  with  electrons  above  50  KeV  and  ions 
above  50  eV. 


Digital  data  (particle  counts)  is  processed  by  use  of  accumula- 
tors. In  all  cases,  the  values  will  be  log  compressed  from  the  16-bit 
accumulators  to  either  10  or  9 bits.  The  program  converts  the  values  on 
tape  from  the  log  compressed  values  back  to  the  particle  counts  before 
outputting  the  results. 

For  the  16  bits,  the  compressor  will  accumulate  the  number  of 

shifts  required  to  reach,  but  not  include,  the  first  "1."  The  number  of 

shifts  required  will  then  be  represented  in  the  9 bit  readout  word  by 

8 5 

the  4 most  significant  bits  (MSB)  of  the  nine  (2  - 2 ) . This  is  called 

the  exponent.  The  compressor  will  then  load  the  next  five  bits  following, 
but  not  including,  the  first  "1"  into  the  five  least  significant  bits 
LSB's  (24  - 2°)  of  the  9 bit  word.  For  10  bit  compression,  six  bits  will 
be  taken  and  the  last,  or  LSB , will  be  loaded  into  another  word  for  readout. 


This  effort  will  be  continued  under  a different  pr  bl  m . : 

Output  will  include  averages,  and  plots  of  the  data,  along  wit 
analysis  of  the  data  tape  supplied. 
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19.  Accuracy  Of  FqF2  Prediction  Techniques 


INITIATOR:  C.  PIKE 

PROBLEM  NO:  4825  PROJECT  NO:  7663 


BACKGROUND 

A copy  of  a program  "ITS  Coefficient  Program,"  was  furnished  for 

modification.  This  program,  along  with  information  obtained  from  the 

ISIS  satellite,  were  part  of  this  analysis  task.  The  objective  of  the 

analysis  was  to  determine  the  accuracy  of  the  ITS  predicted  f F_  and  the 

o 2 

accuracy  of  a second  f^F2  prediction  technique  which  is  based  upon  input 
from  ionograms  obtained  from  the  ISIS  satellite. 

ANALYSIS  AND  RESULTS 

The  program,  "ITS  Coefficient  Program"  calculates  the  following 
coefficients : 


(1) 

f 

o 2 

(2) 

M ( 3000)  F2 

(3) 

h’  . F_ 
min  2 

(4) 

f E 
o 

(5) 

f E (Upper 
os 

Decile) 

(6) 

f E (Medium) 
o s 

(7) 

f E (Lower 
o s 

Decile) 

The  original  purpose  of  the  problem  was  to: 

1)  Determine  exactly  the  logic  of  the  program  to 
facilitate  use  and  modification  of  the  program. 

2)  Document  the  original  program  obtained. 
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3)  Modify  the  program  so  that,  if  desired,  it  will 
produce  only  the  value  of  fQF2  at  a particular 
point,  i.e.,  eliminate  the  worldwide  tables. 

4)  Merge  the  modified  version  of  the  coefficient 
program  to  accept  as  input  the  output  from 
SUA's  program  Ephemeris  Extrapolation  Program 
(LOKANGL) . 

After  all  of  the  above  tasks  had  beer,  completed,  the  modified 
program  was  further  changed  to  accept,  as  input,  either  daily  or  monthly 
sunspot  numbers.  Sunspot  data  had  been  supplied.  To  check  the  program, 
it  has  been  run  to  generate  predictions  for  specific  periods.  These  pre- 
dicted values,  based  both  on  daily  and  monthly  median  sunspot  values,  have 
also  been  prepared  for  future  analysis. 

An  additional  program  was  required  and  written  to  aid  in  the 

analysis  of  determining  the  accuracy  of  the  predicted  fQF2  generated  from 

tne  modified  program  and  the  accuracy  of  an  f F_  prediction  technique  based 

o 2 

on  input  from  the  ionograms.  The  program  reads  in  actual  fQF2  values  and 
values  obtained  from  the  ionograms.  Then  the  program  extracts  and  inter- 
polates when  necessary  from  the  previously  generated  output  the  predicted 
f F froin  the  modified  ITS  Coefficient  Program.  Differences  and  differences 
squared,  are  calculated.  The  values  are  displayed  and  saved  for  future  analysis. 

Periods  of  interest  for  this  analysis  have  been  August,  September, 
November  and  December  1972.  Values  have  been  obtained  for  approximately 
1000  data  points.  These  values,  from  both  prediction  techniques,  will  be 
analyzed  and  the  results  published  in  a report.  (Sample  output  is  shown 
in  Figure  19-1.)  The  analysis  will  consist  of  comparing  the  accuracy  of  the 
two  prediction  techniques.  Preliminary  statistics  will  include  sums  of 
differences  for  various  time  periods  and  rms  values.  Sufficient  data  exists 
to  permit  comparison  among  seasons,  times  of  day,  and  geographic  locations. 

The  modified  ITS  Coefficient  Program,  or  some  version  thereof,  is  being  used 
in  brother  f0F2  analysis  task. 
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Figure  19-1  Sample  f F Output 
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Conversion  Of  Pioneer  Spacecraft  Data 


INITIATOR:  D.  SMART 

PROBLEM  NO:  4930  PROJECT  NO:  2331 


BACKGROUND 

Four  tapes  containing  complete  trajectory  information  of 
Pioneer  6,  7,  8,  and  9 were  received  for  processing.  These  tapes  were 
generated  at  NSSDC  (National  Space  Science  Data  Center)  by  taking  the 
most  accurate  information  from  each  ephemeris  tape  provided  by  JPL  (Jet 
Propulsion  Laboratory)  and  eliminating  time  overlap.  The  data  tapes, 
generated  on  an  IBM  7094,  are  7-track,  800  B.P.I.  Each  logical  record 
contains  89  words,  and  each  block  contains  20  logical  records.  Table 
20-1  is  a listing  of  the  information  that  is  available  on  the  tape.  It 
was  requested  that  these  tapes  be  converted  and  written  on  CDC  SCOPE  tapes 
with  one  logical  record  per  block. 

ANALYSIS  AND  RESULTS 

This  request  was  relatively  simple,  since  a conversion  routine 
had  previously  been  written  to  convert  IBM  7094  tapes.  Section  5 of  this 
report  gives  an  explanation  of  the  conversion.  Changes  were  made  to  pre- 
vious software  by  increasing  the  size  of  the  input  array  and  output  formats. 
The  rest  was  duplicated  from  the  existing  routine. 

There  was  one  problem.  This  was  the  double  precision  time 
word.  The  manner  of  representing  a double  precision  word  was  quite 
different  between  the  two  machines.  It  was  further  noted  that,  because 
of  the  smaller  word  size  (36-bits)  in  the  IBM  7094  computer,  a double 
precision  word  could  almost  be  accommodated  into  the  CDC  60-bit  word. 

This  indicated  a loss  of  12  bits.  However  upon  inspection  of  the  dump, 
it  was  found  that  the  last  6 bits  contained  zeros,  and  precision  was  lost 
in  only  6 bits.  This  was  more  than  was  needed  by  the  requestor. 


This  double  precision  word  was  then  packed  into  one  CDO  6600 
word  and  the  value  output  as  a single  precision  word.  The  four  tapes 
received  were  converted  and  new  tapes  generated  for  the  task  initiator. 
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TABLE  20-1 


PIONEER  SPACECRAFT  DATA 


| WORD  WITHIN 
RECORD 

PARAMETER 

UNIT  OF 
MEASUREMENT 

001 

FIRST  WORD  DOUBLE  PRECISION  TIME 

SEC 

002 

SECOND  WORD  OF  DOUBLE  PRECISION  TIME 

SEC 

003 

JULIAN  DATE 

DAYS 

JULIAN  DATE 

DAYS 

005 

FIRST  WORD  OF  GREG.  CALENDAR  DATE 

YYMMDDDHH 

006 

SECOND  WORD  OF  GREG.  CALENDAR  DATE 

MMSSFFF 

007 

SINGLE  PRECISION  TIME 

SEC 

008 

TIME  FROM  LAUNCH 

SEC 

009 

RADIUS  EARTH  TO  PROBE 

KM 

010 

DECLINATION  OF  PROBE 

DEG 

Oil 

RIGHT  ASCENSION  OF  PROBE 

DEG 

012 

GEOCENTRIC  VELOCITY  VECTOR  MAG 

KM/SEC 

| 013 

INERTIAL  PATH  ANGLE 

DEG 

014 

INERTIAL  AZIMUTH  ANGLE 

DEG 

015 

RADIUS  EARTH  TO  SUN 

KM 

016 

DECLINATION  OF  SUN 

DEG 

017 

RIGHT  ASCENSION  OF  SUN 

DEG 

018 

RADIUS  EARTH  TO  MOON 

KM 

019 

DECLINATION  OF  MOON 

DEG 

020 

RIGHT  ASCENSION  OF  MOON 

DEG 

021 

HELIOCENTRIC  RADIUS  VECTOR  MAGNITUDE 

KM 

022 

CELESTIAL  LATITUDE  OF  PROBE 

KEG 

023 

CELESTIAL  LONGITUDE  OF  PROBE 

DEG 

024 

INERTIAL  VELOCITY 

KM/SEC 

025 

INERTIAL  PT'T’H  ANGLE 

KEG 

026 

CELESTIAL  LATITUDE  OF  EARTH 

DEG 

027 

CELESTIAL  LONGITUDE  OF  EARTH 

DEG 

028 

X COMPONENT  OF  SC  IN  SUN  EARTH  LINE 

KM 

029 

Y COMPONENT  OF  SC  IN  SUN  EARTH  LINE 

KM 

030 

Z COMPONENT  OF  SC  IN  SUN  EARTH  LINE 

KM 

TABLE  20-1  (Continued) 


WORD  WITHIN 
RECORD 

031 

032 

033 

034 

035 

036 

037 

038 

039 

040 

041 

042 
04  3 

044 

045 

046 

047 

048 

049 

050 

051 

052 

053 

054 

055 

056 

057 

058 

059 

060 


PARAMETER 

UNIT  OF 
MEASUREMENT 

RADIUS  SUN  TO  PROBE 

KM 

EARTH  SUN  PROBE  ANGLE 

DEG 

LONGITUDE  OF  SC  IN  SUN  EARTH  LINE  SYS 

DEG 

GEOCENTRIC  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

MOON  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

SUN  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

VENUS  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

MARS  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

SATURN  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

JUPITER  X COORDINATE  OF  SC 

KM 

Y 

KM 

Z 

KM 

GEOCENTRIC  VELOCITY  X COORDINATE  OF  SC 

KM/SEC 

Y 

KM/SEC 

Z 

KM/SEC 

CENTRAL  BODY 

TARGET  BODY 

GEOCENTRIC  LATITUDE 

DEG 

TABLE  20-1  (Continued) 


WORD  WITHIN 
RECORD 

PARAMETER 

— 

UNIT  OF 
MEASUREMENT 

061 

GEOCENTRIC  LONGITUDE 

DEG 

062 

ALTITUDE  ABOVE  EARTH 

KM 

063 

HINGE  ANGLE  OF  EARTH 

DEG 

064 

SWIVEL  ANGLE 

DEG 

065 

HINGE  ANGLE 

DEG 

066 

ALTITUDE  ABOVE  ECLIPTIC  PLANE 

KM 

067 

R(SUN) *COSB 

KM 

068 

R (SUN) *SINB 

KM 

069 

EARTH  PROBE  SUN  ANGLE 

DEG 

070 

EARTH  PROBE  TARGET  ANGLE 

DEG 

071 

EARTH  PROBE  NEAR  LIMB  OF  TARGET  ANGLE 

DEG 

072 

SUN  EARTH  PROBE  ANGLE 

DEG 

073 

SUN  PROBE  NEAR  LIMB  OF  EARTH  ANGLE 

DEG 

074 

SUN  TARGET  PROBE  ANGLE 

DEG 

075 

MOON  EARTH  PROBE  ANGLE 

DEG 

076 

MOON  PROBE  SUN  ANGLE 

DEG 

077 

EARTH  PROBE  MOON  ANGLE 

DEG 

078 

TARGET  PROBE  SUN  ANGLE 

DEG 

079 

CANOPUS  CLOCK  ANGLE  EARTH  CENTER 

DEG 

080 

MOON  CLOCK  ANGLE  EARTH  CENTER 

DEG 

081 

TARGET  CLOCK  ANGLE  EARTH  CENTER 

DEG 

082 

RADIUS  TARGET  TO  SC 

KM 

083 

VELOCITY  SO  WRT  TO  TARGET 

KM/SEC 

084 

ALTITUDE  ABOVE  TARGET 

KM 

085 

ANGULAR  SEMI-DIMETER 

DEG 

086 

EARTH  CLOCK  ANGLE 

DEG 

087 

CANOPUS -PROBE-EARTH  ANGLE 

DEG 

088 

CANOPUS-PROBE-SUN  ANGLE 

DEG 

089 

DAYS  PAST  MIDNITE  FROM  LAUNCH 

DAYS 

95 


21. 


Modification  Of  IEMCAP 


INITIATOR:  P.  ROTHWELL 


PROBLEM  NO:  4915  PROJECT  NO:  7661 


BACKGROUND 

A^GL  PHG  has  received  a copy  of  the  Intrasystem  Electromagnetic 
Compatibility  Analysis  Program  (IEMCAP)  from  RADC . IEMCAP , consisting  of 
two  programs,  models  the  electrical  properties  of  an  aircraft.  The  purpose 
of  this  problem  was  to  get  the  software  supplied  by  RADC  operational  at 
AFGL.  The  software  was  supplied  via  a card  image  tape.  A sample  data  deck 
was  supplied.  It  is  intended  that  the  programs  will  be  used  with  future 
PHG  modeling  efforts. 

ANALYSIS  AND  RESULTS 

The  software  package  obtained  consisted  of  two  programs.  A 
sample  data  deck  was  supplied.  (We  were  not  made  aware  of  the  availability 
of  the  user's  manuals  until  the  task  was  completed.  However  this  did  not 
impede  the  performance  of  the  task.  Rather,  it  gave  us  the  opportunity  to 
investigate  the  logic  of  the  software.)  The  codes  for  the  first  of  the  two 
programs  was  modified  to  allow  for  the  running  at  AFGL.  Calculations  were 
performed  and  the  results  were  checked  to  validate  the  program.  After  vali- 
dating,TART,  the  second  of  the  two  programs,  was  modified  to  allow  for  the 
processing  on  the  CDC  6600.  The  input  cards  required  by  TART  were  determined 
by  investigation  of  the  program.  Further  study  resulted  in  the  matching  of 
the  input  files  required  (and  their  designation)  by  TART  with  output  files 
generated  by  the  first  of  the  two  programs  supplied.  Sample  outputs  were 
generated  for  both  programs  and  results  were  submitted  for  validation.  The 
output  generated  has  been  judged  accurate.  This  software  package  will 
possibly  be  modified  and  incorporated  in  future  modeling  efforts  at 
AFGL/PHG. 
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DESCRIPTION  OF  IEMCAP* 


The  purpose  of  IEMCAP  is  to  determine  whether  signals  from  one  or 
more  emitters  entering  a receptor  port  cause  interference  with  the  required 
operation  of  that  receptor.  Electromagnetic  interference  (EMI)  is  assessed 
in  the  program.  This  is  done  by  computing  "EMI  Margin"  for  each  receptor 
port.  This  EMI  Margin  is  the  ratio  of  power  received  at  a receptor  port  to 
that  receptor's  susceptibility.  The  program  computes  the  margin  in  decibels. 

The  more  positive  the  number  in  decibels  the  greater  is  the  interference 
while  the  more  negative  the  number  in  decibels  the  greater  is  the  compatibility. 

The  program  performs  the  interference  analyses  by  exercising 
various  programmed  formulae  corresponding  the  mathematical  models  of 
emitters,  receptors,  and  the  signal  transfer  mechanisms  between  emitters 
and  receptors.  The  mathematical  models  of  emitters  and  receptors  in  the 
program  correspond  to  the  required  portion  of  their  spectra.  The  pro- 
gram contains  routines  to  compute  the  non-required  spectra  of  emitters 
and  receptors.  These  non-required  spectra  are  initially  based  on  the 
interference  curves  from  military  standards  displaced  in  level  by  the 
user. 


The  following  is  a brief  description  of  the  software  package 
provided.  The  first  step  in  using  IEMCAP  is  to  assemble  the  appropriate 
data  for  the  system  to  be  analyzed  on  punched  cards.  This  data  is  then 
fed  into  the  Input  Decode  and  Initial  Processing  Routine  (DIPR)  section 
of  IEMCAP. 


* For  additional  information  on  the  basic  mathematical  functions 
required,  see  Intrasystem  Electromagnetic  Compatibility  Analysis  Program, 
Volumes  I and  II,  McDonnell  Aircraft  Corporation,  December  1974,  RADC-TC-34L. 


The  first  section  of  IDIPR  decodes  the  punched  cards  and  checks 
for  errors.  If  an  error  is  detected  a diagnostic  message  is  printed,  the 
card  is  deleted, and  the  program  continues  processing  the  rest  of  the  data. 

The  program  normally  stops  after  all  data  has  been  read  if  there  were  errors. 

If  there  were  no  input  card  errors,  the  Initial  Processing  Routine 
(IPR)  section  of  IDIPR  is  executed.  For  the  first  run  of  a system,  the 
spectrum  math  models  are  accessed  for  each  port.  These  models  use  the 
user-supplied  spectrum  parameters  to  generate  the  spectra.  The  processed 
user  and  spectrum  data  is  written  on  a disk  file  called  the  Intrasystem 
File  or  Intrasystem  Signature  File  (ISF) . This  ISF  can  be  used  as  input 
for  subsequent  runs,  either  as  is  or  modified  by  additional  card  inputs. 

Also,  a new  ISF  can  be  generated  containing  the  modifications.  Thus,  an 
updated  ISF  can  be  maintained  for  changes  in  system  design.  The  data  is 
also  written  on  a number  of  working  files  for  use  by  the  Task  Analysis 
Routine  (TART)  program  of  IEMCAP,  and  a printed  report  of  the  data  is 
generated. 

The  TART  section  of  IEMCAP  performs  four  basic  analysis  tasks. 

These  tasks  are  summarized  below: 


• Specification  Generation  - Adjusts  the  initial 
non-required  emission  and  susceptibility  spectra 
for  system  compatibility.  The  user-specified 
adjustment  limit  prevents  too  stringent  adjust- 
ments. A summary  of  interference  situations  not 
controllable  by  EMC  specif ications  is  displayed. 

The  adjusted  spectra  are  the  maximum  emission  and 
and  minimum  susceptibility  specifications  for  use 
in  EMC  tests.  Analysis  results  are  stored  on  the 
Baseline  Transfer  File  (BTF ) for  subsequent  runs. 

• Baseline  System  EMC  Survey  - Surveys  the  system  for 
interference.  If  maximum  EMI  margin  over  the  fre- 
quency range  for  a coupled  emitter-receptor  port 
pair  exceeds  the  user-specified  printout  limit,  a 
summary  of  the  interference  is  printed.  Total 
received  signal  into  each  receptor  from  emitters 

is  printed.  The  results  are  stored  on  the  BTF 
for  subsequent  runs. 
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• Trade-off  Analysis  - Compares  the  interference 
for  a modified  system  with  values  stored  on  the 
BTF  from  previous  run.  This  is  to  determine 
the  effect  on  interference  of  antenna  changes, 
filter  changes,  spectrum  parameter  changes,  wire 
changes,  or  other  changes. 

• Specification  Waiver  Analysis  - Shifts  specified 
portions  of  specific  port  spectra  and  compares  the 
resulting  interference  to  that  stored  on  the  BTF. 
The  effect  of  granting  waivers  for  specific  ports 
can  be  assessed. 


TART  is  composed  of  two  basic  routines.  The  Specifications  Genera- 
tion Routine  (SGR)  which  performs  the  first  task  above,  and  the  Comparative 
EMI  Analysis  Routine  (CEAR)  which  performs  the  other  three.  These  access 
the  coupling  math  model  routines  to  compute  the  transfer  ratios  between 
emitter  and  receptor  ports.  The  two  parts  of  IEMCAP  are  executed  separately 
with  data  files  used  for  intermediate  storage  between  parts.  Central  Pro- 
cessing Unit  (CPU)  core  memory  to  load  and  execute  each  part  of  IEMCAP  on 
CDC  6600  at  AFGL  are  as  follows: 

IDIPR  - 235K  words  (octal) 

TART  - 225K  words  (octal) 

The  program  uses  4 permanent  and  10  working  files,  in  addition  to 
normal  card  input  and  printed  output.  The  permanent  files  are  presently 
disk  files  at  AFGL.  For  the  working  files,  disk  is  used  because  of  the 
number  of  files  and  large  number  of  accesses  per  run.  All  files  are 
sequential . 

The  amount  of  file  space  needed  depends  on  the  size  of  the  system 
being  analyzed.  Except  for  the  BTF,  the  size  depends  primarily  on  the 
number  of  ports.  The  size  of  the  BTF  depends  on  the  number  of  coupled 
port  pairs  since  it  stores  the  analysis  results. 

The  execution  time  also  depends  on  the  system  size  and  analysis 
task  desired.  TART  run  time  primarily  depends  on  the  number  of  coupled 
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port  pairs.  This  potentially  increases  as  the  square  of  the  number  of 
ports.  However,  each  emitter  port  is  not  generally  coupled  to  each 
receptor  port  so  the  actual  time  will  be  less.  Specification  generation 
requires  three  passes  through  the  emitters  per  receptor  and  two  passes 
through  the  receptors  per  run  and  runs  longer  than  the  other  tasks. 

DATA  STORAGE  FILES 

As  discussed  above,  IEMCAP  uses  a number  of  data  files.  Each 
physical  file  and  device  is  assigned  a logical  unit  number  which  is 
stored  in  each  section  of  the  program  as  a mnemonic  variable  name.  For 
example,  on  the  CDC  6600  the  card  input  is  designated  as  logical  unit 
number  5.  In  IDIPR  and  TART,  the  variable  INN  is  set  to  5,  and  all  card 
input  read  statements  reference  INN.  The  logical  unit  numbers  are 
assigned  to  files  in  the  job  control  cards.  Also,  a number  of  files  are 
assigned  to  the  same  logical  unit.  This  allows  multiple  usage  of  the  same 
physical  file  space. 

The  files  are  categorized  as  permanent,  working  and  scratch. 
Permanent  files  are  used  to  store  data  and  analysis  results  for  use  in 
subsequent  runs.  Working  or  intermediate  files  provide  temporary  storage 
of  the  data  for  efficient  use  by  the  various  routines.  They  also  provide 
intermediate  data  storage  between  IDIPR  and  TART.  Scratch  files  are  used 
for  temporary  storage  within  IDIPR  and  TART. 

The  discussion  in  the  paragraphs  immediately  following  is  a 
discussion  of  the  system  mathematical  model  employed  in  IEMCAP. 

SYSTEM  MATHEMATICAL  BASIS 


The  basis  for  all  calculations  performed  by  each  of  the  functions 
of  IEMCAP  is  the  linear  relationship  for  power  coupled  from  an  emitter, 
through  a transfer  medium,  and  received  by  a receptor.  The  general 
communication  theory  equation  relating  power  spectral  density  at  the 
detector  of  a receptor  to  power  spectral  dersity  present  at  an  emitter's 
output  port  is  expressed  as  follows: 
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(1) 


PSD 


where 


PSD 


Vf> 

6.  (f> 

1 


= r>  (f ) T.  .(f)  B.  (f) 
S]  13  l 


= output  power  spectral  density  (in  watts/Hz) 
received  by  receptor  i (at  its  detector)  from 
emitter  j, 

= output  power  spectral  density  (in  watts/Hz) 
at  the  terminals  of  source  j (including  cw 
power  as  delta  functions) , 

= power  transfer  function  of  the  coupling  medium 
between  source  j and  receptor  i, 

= receptor  response  function  relating  power  at 
the  detector  to  power  at  the  input  terminals. 


The  broadband  and  narrowband  components  of  the  received  power 
spectral  density  may  be  computed  separately  by  considering  the  broad- 
band and  narrowband  components  of  the  output  power  spectral  density  of 
the  source  j.  Thus,  the  power  spectral  density  of  source  j is  expressed 

33 ; 


Vf)  = 


nSB3(f) 


(2) 


where 


n • (f) 

SB  j 


the  broadband  power  spectral  density  of  source  3 


and 


nSNj(f) 


the  narrowband  power  spectral  density  of  source  j 
(composed  of  delta  functions  at  the  frequencies 
of  the  cw  signals.) 
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Spectral  Analysis 


INITIATOR:  T.  Conley 


PROBLEM  NO.:  4865  PROJECT  NO.:  7670 


BACKGROUND 


ASEC  is  currently  generating  software  to: 

1.  Modify  a real  or  synthetic  digital  input  spectrum 

by  an  appropriate  atmospheric  transmittance  function; 

2.  Convolve  it  with  an  instrument  function; 

3.  Process  and  plot  the  resultant  spectrum.  The  digital 
processing  scheme  provides  the  framework  for  com- 
parison of  synthetic  spectra,  infrared  backgrounds, 
integrated  measurements,  and  synthetic  spectra  genera- 
tion models. 


The  following  describes  what  has  been  accomplished  to  date. 


ANALYSIS  AND  RESULTS 

Convolution  of  two  functions  is  an  important  operation  in  many 
scientific  fields.  In  particular,  this  analysis  of  measured  infrared  spec- 
tra convolves  the  synthetic  spectrum  with  an  instrumer'  function. 

The  convolution  integral,  which  is  the  basis  of  the  many  con- 
volution methods,  is  given  by: 

00 

Y(t)  = J X(T)  h (t  - T)  dT 
= X(t)  *h(t) 
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That  is,  the  convolution  of  the  functions  X(t)  and  H(t)  yields  the  func- 
tion Y (t)  . 

The  procedure  for  convolving  by  the  integral  method  is  as  follows 

1.  Take  the  mirror  image  of  h(T)  about  the  ordinate  axis. 

2.  Shift  h(-T)  by  the  amount  t. 

3.  Multiply  the  shifted  function  h (t  - r)  by  X(x). 

4.  Area  under  the  product  of  h(t  - T)  and  X(T)  is  the 
value  of  the  convolution  at  time  t. 

The  general  convolution  limits  are  -«  to  +<».  It  is  desired  to 
find  the  lower  and  upper  limits  so  that  they  will  contain  limits  of  both 
functions.  To  do  so,  we  choose  the  smaller  of  the  two  lower  limits  and 
the  larger  of  the  two  upper  limits.  Tht se  two  limits  will  give  us  the 
common  least  upper  bound  and  greatest  lower  bound  for  both  functions. 

Applications  of  the  Fast  Fourier  Transforms  (FFT)  in  digital 
filtering,  power  spectrum  analysis,  simulation,  etc.,  are  based  on  a 
specific  implementation  of  the  convolution  integral  which  employs  the 
FFT  as  an  approximation  to  the  continuous  Fourier  transform.  The  FFT  is 
simply  a procedure  for  rapidly  computing  the  discrete  Fourier  transform. 
It  is  advantageous  to  use  the  FFT  in  computing  the  convolution  due  to  its 
tremendous  increase  in  computational  speed. 

The  FFT  convolution  of  finite  duration  for  the  discrete  case 
may  be  defined  as  follows: 

N-l 

Y(K)  = XU)  h(K  " L) 
i=0 

where  both  X (K)  and  h(K)  are  periodic  functions  with  period  N.  The 
discrete  convolution,  if  correctly  performed,  will  produce  a replica 
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of  the  continuous  convolution  provided  that  both  the  functions  X(t) 

and  h(t)  are  of  finite  duration.  Evaluation  of  the  N samples  of  the 

convolution  using  the  above  equation  requires  a computation  time 
2 

proportional  to  N (the  number  of  multiplications). 

Special  care  must  be  taken  if  the  two  functions  have  different 
periods  or  if  at  least  one  of  the  functions  has  unevenly  spaced 
samples.  Impulse  functions  are  of  the  unevenly  spaced  type.  It  is 
important  to  make  sure  that  the  function  being  shifted  can  be  evaluated 
for  all  values  of  the  impulse  function. 

On  Contract  F19628-76-C-0203,  ASEC  developed  a special  routine 

to  convolve  such  functions.  The  problem  was  to  convolve  the  intensity 

/ 2 

function  with  an  instrument  function  approximated  by  ((sin  X)/X)  . 

In  theoretical  discussions  the  wave  number  v = v^/c  = 1/A  is 
generally  used.  (A  is  wave  length  of  the  spectral  lines  in  the  infra- 
red; c is  speed  of  light  in  a vacuum;  = c/A  is  frequency).  Wave  number 
is  proportional  to  the  frequency  and  is  measured  in  cm  1 (number  of  waves 
per  cm) . The  frequency  and  wave  number  depend  on  the  energy  E according 
to  Planck's  relation  E = hv^  = hoc,  where  h is  Planck's  constant.  According  to 
the  Maxwell-Boltzmann  distribution  law  the  number  of  molecules  that  have 

a classical  vibrational  energy  between  some  E and  E + AE  is  proportional 
-E/kt 

to  e dE  where  k is  Boltzmann’s  constant  and  T is  the  absolute 

temperature.  In  quantum  theory,  only  discrete  values  are  possible  for 

vibrational  energy.  The  number  of  molecules  in  each  vibrational  state 

-E/kT 

is  again  proportional  to  the  Boltzmann  factor  e . The  zero-point 

energy  can  be  left  out. 

The  wave  number  v of  a spectral  line  can  be  represented  as 
the  difference  between  two  quantities  called  terms.  Terms  are  negative 
energy  values  divided  by  hC.  They  can  be  obtained  empirically  for  the 
observed  spectrum.  The  first  term  is  always  given  by  the  wave  number  of 
the  limit  of  the  series  of  spectral  lines  in  question.  Often  the  word 
term  is  used  interchangeably  with  energy  state  or  quantum  state. 


In  order  to  determine  the  total  number  of  molecules  N,  one  must 
consider  that  N is  proportional  to  the  sum  of  the  Boltzmann  factor  over 
all  states,  i.e.,  the  state  sum  of  partition  function. 


The  number  of  molecules  in  the  state  v is 

N = JL  e *Go(v)hc/kt 
v Qv 


where 


+ e~G0 (l)hc/kt+  c-G0(2)hc/kt  + 


and 


-0o (v) hc/kt 


-E/kt 

e 


for  state  v 


The  probabilities  of  transition  under  the  influence  of  radiation 
are  determined  by  the  eigenfunctions  of  the  states  involved.  The  in- 
tensities of  the  emitted  or  absorbed  spectral  lines  can  be  determined. 
When  the  eigenfunctions  are  known  one  can  calculate  whether  or  not  two 
states  can  combine  with  each  other.  Transitions  that  cannot  occur  are 
considered  forbidden  transitions. 


The  intensity  of  a spectral  line  in  emission  is  defined  as  the 

energy  emitted  by  the  source  per  second.  Given  atoms  in  the  initial 

state  and  A fractions  of  atoms  in  the  initial  state  carrying  out  the 
nm 

transition  to  m per  second  then  I = N hcv  A 

em  n nm  nm 


hcv  is  the  energy  of  each  light  quantum  of  wave 

numBer  v emitted  in  transition, 
nm 


A is  the  Einstein  transition  probability  of  spontaneous  emission. 

The  intensity  of  absorption  is  given  by  in[*  = p N 6 Axhcv  where  N 

abs  nm  m mn  nm  m 

is  the  number  of  atoms  per  cm  in  the  lower  state  m and  B^  is  the  Einstein 
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transition  probability  of  absorption.  p N 3 represents  the  number  of 
^ nn  m mn 

transitions  per  cm  per  second  produced  by  incident  radiation,  AX  is  the  thick- 
ness of  the  layer.  For  this  task  we  were  primarily  interested  in  the  absorption 
intensi ties. 


In  the  software  the  partition  function  values  were  computed  using 
the  given  equation, 

15.5 

(2J  + 1)  EXP 


Q (T) 
*v 


E 

J= . 5 


=§  VJ> 


where 


Ev<J)  = rotational  energy 


H = 3.3356  x 10  11  cm  1 sec  (Plank  constant) 


C = 3 x 1010  cm/sec  (speed  of  light) 


K = 6.95  x 10  1 cm  1 °K  1 (Boltzmann  constant) 


T = 220°K  (Temperature) 


E (J)  = J ( J + 1)BV,  where  Bv  is  an  array  of  constant  values 
corresponding  to  molecules  of  interest 


J = total  angular  momentum  (rotational  level) 


Second,  intensities  were  computed  by: 

(2J 


N ( J , V)  = N 


v Q. 


-7=7^-  EXP  ^ E (J) 
(T)  KJ  v 

f u J 


where  is  an  array  of  population  values. 
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Figure  22-1  Plot  of  Normalized  Intensity  Versus  Wave  Number 


A sample  of  normalized  intensity  versus  wave  number  is  a i ,en  in 

Figure  22-1.  The  intensity  values  were  then  convolved  with  tie  instrument 
2 

function  (sin  X/X)  in  the  following  manner: 

1.  Intensity  function  was  shifted  by  the  resolution 
value  and  rotated  about  the  resolution  value. 

2.  The  new  function  was  shifted  to  the  tight  by  dX 
and  the  instrument  function  (sin  X/X)^  was 
evaluated  at  the  intensity  valuer,  in  the  range 
of  the  resolution. 

3.  These  two  functions  were  multiplied  and  summed 
over  the  product  values  to  give  one  convolution 
value . 

The  above  steps  were  repeated  until  the  left-most  value  of  intensity 
function  was  out  of  the  range  of  the  instrument  function. 

At  the  end,  the  convolution  values  versus  the  wave  number  were 
plotted.  The  plotted  values  as  well  as  the  digital  values  are  used  in 
comparison  of  synthetic  spectra,  infrared  backgrounds  and  integrated 
measurements.  Present  software  provides  a basis  for  the  presentation 
of  the  synthetic  spectra  generation  models.  The  original  software  has, 
in  the  past,  and  presumably  will  in  the  future,  be  adopted,  modified 
and/or  expanded  to  perform  continued  analyses. 
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23.  Signal  Statistics  Of  Scintillation 


INITIATOR:  H.  WHITNEY 

PROBLEM  NO:  4893  PROJECT  NO:  4643 


BACKGROUND 

F-layer  irregu’  a rities  in  the  ionosphere  can  cause  variations  in 
radio  signals  which  traverse  a disturbed  region  of  the  ionosphere.  These 
signal  variations  are  known  as  ionosphere  scintillations.  Performance  of 
satellite  communications  at  UHF  will  be  affected  when  scintillations  exceed 
the  fade  margin.  The  effect  is  most  noticeable  when  propagation  paths  cross 
through  the  auroral  and  equatorial  ionosphere. 

Ionosphere  scintillation  data  are  being  analyzed  by  ASEC  on  the 
current  contract.  Software  was  written  and  digitized  data  tapes  were 
processed.  The  results  were  used  for  comparison  with  those  of  earlier 
studies.  Hie  digitized  sampling  rate  of  the  data  was  varied  from  six 
samples/second  to  36  samples/second  to  determine  the  effect  on  the  analysis. 

ANALYSIS  AND  RESULTS 

Data  from  Peru  is  being  used  for  this  study.  The  signal  amplitude 
from  satellite  beacons  was  recorded  on  FM  analog  tape  and  then  digitized. 
Calibration  of  equipment  was  accomplished  and  recorded  prior  to  and  sub- 
sequent to  acquisition  of  data.  Computer  programs  were  written  to  convert 
the  digitized  values  into  readable  form  by  taking  the  data  values  from  the 
original  tape  and  writing  them  out  on  another  tape  in  CDC  6600  format  words. 
The  original  tapes  contain  both  calibration  files  and  data  files. 

After  the  data  was  written  onto  CDC  tapes,  the  calibration  files 
are  used  to  convert  the  digitized  data  into  signal  amplitudes  in  dB.  All 
dB  numbers  are  negative,  with  0 as  the  maximum.  The  dB  values  are  written 
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onto  magnetic  tapes  for  permanent  storage,  and  onto  permfiles  for  use  with 
the  analysis  programs. 

Data  was  originally  analyzed  in  15-minute  intervals.  In  addition, 
the  cumulative  distribution  was  calculated  on  1.5-minute  intervals,  as  well 
as  for  15-minute  intervals,  to  determine  if  a shorter  interval  size  is  more 
representative  of  broad  changes  in  the  data.  Plots  of  the  data  were  obtained 
for  each  1.5-minute  interval  (Figure  23-1)  and  for  cumulative  amplitude 
distribution  (Figure  23-2) . 

Software  used  for  the  analysis  represents  a combination  of 
both  totally  new  routines  and  adaptation  of  existing  ones.  A program 
was  written  to  calculate  a number  of  different  statistical  relation- 
ships for  the  given  scintillation  data.  First,  the  data  in  dB  levels, 
is  plotted  on  a graph,  giving  dB  versus  time  in  seconds,  for  each 
channel . 


Second,  the  autocorrelations  of  data  in  power  are  found  for  each 
channel,  for  0 - to  16-second  time  lags,  at  intervals  of  At  = 1/6  second. 

The  dB  levels  (DJ  are  converted  to  power  (WJ  by  using  the  conversion  formula 


W. 

l 


10 


(10  Di) 


These  autocorrelations  are  plotted  versus  time  lag.  The 
subroutine  FTAUTO  is  used  to  find  the  autocorrelations,  with  the 
formula : 


AC  . 
3 


- T 

N Lu 


W) 


W) 


for  the  Jth  autocorrelation. 


(W  = mean  of  time  series  W,  o 


2 


variance  of  W) 
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Figure  23-1  Sample  Peru  Data 


Figure  23-2  Cumulative  Amplitude  Distribution 


Figure23-3  is  a sample  of  the  autocorrelation  coefficients 

plotted  with  respect  to  the  correlation  interval  for  the  first  16 

seconds.  Autocorrelation  coefficients  greater  than  or  equal  to  0.0  are 

printed  together  with  their  corresponding  lag  times.  The  program  for 

calculating  the  coefficients  first  evaluates  the  mean  of  the  data  value 
— 2 

X and  the  variance  a before  computing  the  coefficients  AC. 

The  autocorrelation  is  a method  of  characteristizing  the  rate 
of  scintillation  fading.  The  width  of  the  autocorrelation  is  inversely 
related  to  the  variance  of  the  data.  Large  variations  at  a fast  rate 
w 11  produce  a short  width  in  the  autocorrelation  function. 


A third  routine  for  this  analysis  calculates  the  crosscorrelation 
between  the  data  in  power  of  the  two  channels  for  time  lags  of  0 to  1 
minute,  with  channel  1 leading  for  the  positive  time  lags  and  channel  2 
leading  for  the  negative  time  lags  The  subroutine  uses  the  formula: 


N-J 


CC 


N / J 


'X.  - X)  (Y  - Y) 

1 + 3 


1=1 


J~2  2 

a a 


for  the  Jth  crosscorrelation.  (X,  Y are  the  means  of  the  two  time  series; 

2 2 

0 , 0 are  the  variances.)  The  crosscorrelations  are  calculated  by  calling 
x y 

the  crosscorrelation  routine  twice,  first  with  channel  1 in  the  data  array, 
the  second  time  with  channel  2 in  the  array.  The  first  value  in  the  out- 
put array,  AC,  is  the  crosscorrelation  at  time  lag  zero;  this  value  is 
plotted  in  the  center  of  the  graph.  The  other  time  lags  are  plotted  to 
the  right  (positive  time  lags,  channel  1 leading)  and  to  the  left  (negative 
time  lags,  channel  2 leading). 


The  cross  correlation  function  is  shown  in  Figure  23-4  for  time 
lags  of  up  to  1.5  minutes  with  first  one  time  series  leading  the  other, 
and  vice  versa.  The  maximum  crosscorrelation  coefficient  is  printed  out. 
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The  fourth  routine  in  the  analysis  software  computer  and  plots  the 
probability  density  function  (PDF)  and  the  cumulative  amplitude  distribu- 
tion function  for  the  data,  in  dB . The  PDF  is  plotted  directly  from  the 
data  by  findinq  the  percentaqe  of  points  at  each  dB  level.  The  Rayleigh 
PDF  is  plotted  on  the  same  graph  for  comparison. 

Tlie  CDF  is  computed  from  the  data  by  findinq  the  cumulative  sum 
of  the  PDF's  for  each  dB . It  is  then  plotted  on  a graph  together  with  the 
Rayleigh  CDF  for  comparison.  The  graphs  are  done  on  a probability  scale 
for  the  X-axis,  and  a linear  scale  in  dB  only  from  -16  to  +8  is  displayed 
for  the  Y-axis.  A probability  of  50%  corresponds  to  dB  = 0. 

Figure  23-1  is  representative  of  strong  scintillations;  and,  as 
shown  by  the  CDF  (Figure  23-2),  its  distribution  closely  resembles  a Ray- 
leigh Distribution.  The  scintillation  index  S , defined  by  Briggs  and 

1 4 
Parkins  , is  calculated  for  each  15-minute  sample.  is  printed  on  all 

plots.  This  value  is  inversely  related  to  the  Nakagami-m  value  (m  = 1/S^) 

and  is  the  most  widely  used  index.  It  is  related  to  the  power  P and  defined  by 

-2  -2 
c2  - P - P 

4 _ 2 

P 


A Chi-Square  goodness  of  fit  test  was  applied  to  the  experimental 

data  to  determine  if  the  distribution  of  the  data  may  be  better  categorized 

as  Log-Normal  or  as  Nakagami-m  (Rayleigh  Distribution  being  a special  class 

. 2 

of  the  Nakagami-m  distribution,  when  m = 1) . The  program  utilizing  the 
Chi-Square  algorithm  tests  how  well  the  real  data,  after  the  conversion 
to  power,  fits  a Rayleigh  Distribution,  Normal  Distribution,  or  a Nakagami-m 
distribution. 

^Briggs,  G.H.,  and  Parkin,  I. A.  (1963),  On  the  Variation  of  Radio 
Star  and  Satellite  Scintillations  with  Zenith  Angle,  J.  Atms.  Terr.  Phys., 

25 : PP • 330- 366 . 

2 

Nakagami,  M (1960),  Statistical  Methods  in  Radio  Wave  Propagation, 
edited  by  W.C.  Hoffman,  PP.  3-36,  Pergamon,  New  York. 
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The  real  data  is  compared  with  a theoretical  distribution  func- 
tion. The  probability  that  the  data  represents  that  particular  distri- 
bution is  calculated.  The  theoretical  distribution  function  is  supplied 
by  the  user.  Equi-probable  cells  are  set  up,  and  the  number  of  data  points 
which  fall  into  each  cell  are  counted.  The  expected  value  of  each  cell  is 
compared  with  the  actual  value  of  each  cell,  and  the  differences  are  saved. 
The  sum  of  the  differences  is  used  to  calculate  the  Chi-Square  statistic, 
which  is  then  used  to  find  the  probability  that  the  real  data  fits  the 
theoretical  distribution.  This  particular  application  uses  25  equi-probable 
cells,  with  22  degrees  of  freedom.  The  steps  in  the  Chi-Square  test  are 
as  follows: 


(1)  Determine  k,  the  number  of  cells 

(2)  Determine  the  hypothesized  distribution,  F(X) 

(3)  Find  the  values  of  X.  such  that  if  the  true 
distribution  is  F(X)  then 


F (xi  ♦ 1>  - F (V  = k 


where 


X = -a,  X = °° 
o k 

X.  determines  the  k cells 

l 

(4)  Calculate  the  number  of  observation  in  the 
cells  called  , for  i = 1,  ...,  k 

(5)  Compute  the  Chi-Square  statistic 

X2  = Z(N.  - r,/k)2/(n/k) 

-2 

(6)  If  F (X)  is  the  correct  distribution  X has 

a Chi-Square  distribution  with  k - 1 degrees 
of  freedom. 

When  F (X)  has  estimates  rather  than  the  actual 
parameter  the  distribution  of  X2  is  Chi-Square 
distributed  with  K - m - 1 degrees  of  freedom 
when  F (X)  is  the  true  distribution. 
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The  results  of  the  Chi-Square  goodness  of  fit  test  are  printed 
out,  along  with  the  mean,  the  median,  the  Nakagami-m,  the  moments,  S^, 
B^,  B.,,  , and  the  values  of  the  CDF  function. 

£ ,Yl  ' " 

M = = moments 

k n 

— = mean  of  data 

y 

y ^ = data  points 
n = number  of  data  points 


/b^  = skewness 
M, 


' B 


1 


<M,) 


3/2 


B = Kurtosis 

M4 

B2  2 
<m2) 


A fifth  program  calculates  the  power  spectrum,  which  is  a 
fast  Fourier  transform  of  the  power  levels  of  the  input  data.  The  output 
is  plotted  on  a semi-log  graph,  on  which  the  Y-axis  is  linear  in  dB  values, 
and  the  X-axis  is  logarithmic  in  frequency  values.  The  ith  frequency  is 
calculated  using  the  formula 

i 

Fi  = 2N  (At) 

where  N = number  of  values  in  the  power  series,  and  At  = sample  rate  in 

seconds.  Optionally,  the  data  series  may  be  split  into  N shorter  segments, 
for  each  of  which  the  spectral  computations  are  made  separately.  The  N 
spectra  are  then  averaged.  For  most  of  the  analysis  N = 1,  so  that  there  is 
only  a single  segment  to  accomodate  an  entire  input  series.  The  number  of 
terms  In  that  series  must  be  a power  of  2.  If  this  is  not  the  case,  the 
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excess  is  zero-filled,  up  to  a maximum,  which  in  our  present  case  is 
8192  words.  The  equation  for  the  power  spectrum  is 


PS.  = 
1 


At 


(X2 

1 


+Xi  ♦!> 


where  PS^  = the  ith  value  of  the  power  spectrum,  Ax  = sample  rate,  N = the 
number  of  points,  = power  series  value.  A ten-point  running  average  is 
taken  on  the  spectral  values,  after  which  they  are  converted  back  to  dB 
values  for  plotting. 


Figure  23-5  is  a sample  of  a power  spectrum.  Various  data  windows 
and  segmentation  may  be  applied  to  smooth  the  incoming  signal  in  order  to 
obtain  more  information  for  future  detailed  analysis.  However,  here  neither 
course  was  taken.  We  chose  instead  to  take  a ten-point  running  mean  of 
the  power  spectra  values  after  they  were  computed. 


The  power  spectrum  is  a useful  way  of  characterizing  the  rate  of 
scintillation  fading.  For  strong  scintillations,  the  power  spectrum  is 
relatively  flat  at  low  frequencies  and  exhibits  a roll-off  at  higher 
frequencies.  The  cut-off  frequency  denotes  the  bandwidth  of  the  process. 

This  parameter  is  related  to  the  inverse  of  the  width  of  the  autocorrelation. 
Autocorrelations  which  have  relatively  large  widths  will  reflect  in  the 
power  spectrum  large  bandwidth  with  no  cut-off  frequency  or  roll-off.  The 
power  spectrum  will  then  appear  to  contain  only  one  measurably  slope. 


There  is  a control  card, read  in  at  the  beginning  of  the  program, 
which  controls  the  sections  of  the  program.  The  first  ten  columns  in  the 
card  are  assigned  to  ten  control  words  • A 0 in  a column  skips  the  corres- 
ponding part  of  the  program ;a  non-zero  directs  action  to  be  taken.  The 
items  controlled  are  as  follows: 

Column  1:  plot  data,  channel  1. 

2:  plot  data,  channel  2. 

calculate  and  plot  autocorrelations,  channel  1. 
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3: 


4:  calculate  and  plot  autocorrelations,  channel  2. 


5:  calculate  and  plot  crosscorrelations. 

6:  calculate  and  plot  power  spectrum,  channel  1. 

7:  calculate  and  plot  power  spectrum,  channel  2. 

8:  calculate  PDF  and  CDF,  channel  1 and  print, 

along  with  results  of  Chi-Square  test. 

9:  calculate  PDF  and  CDF,  channel  2 and  print, 

along  with  results  of  Chi-Square  test. 

10:  plot  PDF  and  CDF. 

Presently  a separate  program  is  being  used  to  calculate  eight 
different  values  which  are  used  to  analyze  the  data.  The  values  can  be 
computed  for  either  1-1/2  minute  segments  or  15  minute  segments.  They 
are  printed  out  at  the  end  of  the  program,  first  for  channel  1 and  then 
for  channel  2.  The  eight  values  are:  the  median  of  the  data  (in  dB); 
the  mean  of  the  power;  S^ ; Nakagami-m;  the  time  for  which  the  auto- 
correlation first  is  equal  to  or  less  than  0.5;the  maximum  crosscorrelation; 
the  time  for  which  the  crosscorrelation  is  at  a maximum;  and  the  velocity, 
i.e.,  366  meters/time  of  maximum  crosscorrelation.  The  autocorrelation 
and  crosscorrelation  are  calculated  from  power  levels.  This  program  will 
be  incorporated  into  the  main  program. 

From  the  observed  signal  levels,  the  amplitude  and  rate  character- 
istics of  intense  scintillations  were  analyzed.  It  can  be  seen  from  the 
distribution  that  the  data  closely  resemble  a Rayleigh  distribution  which 
is  a special  case  of  the  Nakagami-m  distribution  for  m = 1.  The  auto- 
correlation and  power  spectrum  define  the  fading  rate  and  give  a basis 
for  evaluation  of  different  data.  The  crosscorrelation  gives  a means  of 
evaluating  space  diversity  at  similar  frequencies.  Confidence  bands  are 
included  to  give  credence  to  the  reliability  of  the  results.  Also  to  be 
determined  are  (a)  the  effect  a different  digitized  sampling  rate  would 
have  on  the  results,  and  (b)  the  optimal  time  interval  of  analysis.  Addi- 
tional data  exists  at  present  and  more  data  is  being  obtained.  To  obtain 
this  analysis,  all  data  will  be  processed;  and  tables  of  the  output  will  be 
assembled  along  with  plots  for  comparison  of  data. 


121 


24.  Error  Analysis  of  Turbulence  Parameters 


INITIATOR:  E.  MURPHY 

PROBLEM  NO:  4889  PROJECT  NO:  6687 


BACKGROUND 

Seven  years  of  rocket  grenade  data  (winds  and  temperature) , 
representing  153  measurements  from  four  sites  at  different  latitudes, 
have  been  collected.  This  data  has  been  analyzed  statistically  to  pro- 
vide probability  of  occurrence  models  for  turbulence  and  average  value 
of  turbulence  parameters  as  a function  of  latitude,  altitude,  and  time 
of  year.  The  experimental  errors  are  available  for  the  wind  velocity 
and  temperature  data  used.  Equations  or  expressions  which  use  temperature 
and  component  winds  and  their  derivatives  were  supplied  along  with  the 
data  and  the  experimental  errors. 


ANALYSIS  AND  RESULTS 


The  following  equation  is  used  to  determine  the  Richardson  num- 
ber (R^) . The  Richardson  number  is  the  basis  for  the  criterion  for  the 
presence  or  absence  of  turbulence. 

*i  - * H 


where  g = the  gravitational  value  at  the  latitude  - 
assumed  to  have  negligible  error 

T = dry  adiabatic  lapse  rate  which  is  a generally  accepted 
constant  (9.8°k/km) 


3t  ,3v. 2 3Vx. 2 3Vy.  2 

The  remaining  parameters  (T,  ^ and  (^)  = (^— ) + (gp-)  ) are  obtained 

from  the  grenade  data  using  a spline  fit  which  provides  the  temperature, 

3t.  3vy  3V„ 

temperature  gradient  (gp  , and  the  wind  component  gradients  and  -gp ) 

at  the  desired  altitudes. 


r 
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When  R^  < ...5,  turbulence  is  considered  present  and  turbulence 
parameters  are  calculated  using  the  following  relations. 


where 


*'V1/2  * b 


a = 


+0.15  for  R.  > 0 

l 

|-0.15  for  Ri  < 0 


b = 0.08 
V = wind  speed 
N 


+ r)/Ti  1/2 


2 (W' ) N 

e = — — - — — rate  of  dissipation  of  KE 
(W')2  = turbulent  intensity 


K = C 


W 


N 


turbulent  diffusivity 


These  parameters  are  calculated  and  averaged.  They  assume  a value  of  zero 
when  there  is  no  occurrence  of  turbulence. 


The  first  concern  was  to  provide  meaningful  expressions  for  the 
errors  in  the  turbulence  parameters  (rate  of  dissipation  of  kinetic  energy, 
turbulent  diffusivity,  and  turbulent  intensity) , and  to  display  these  re- 
sults using  the  CALCOMP  plotter.  Standard  deviations  were  calculated  for 
the  turbulence  parameters,  the  Richardson  numbers,  and  also  for  temperature 
and  component  winds.  In  addition,  software  was  written  to  introduce 
maximum  error  into  the  calculations.  This  was  done  by  applying  the  errors 
to  the  wind  and  the  temperature  data  in  such  a way  as  to  produce  the 
greatest  errors  in  a spline  interpolation,  and  therefore  in  the  Richardson 
number.  Twenty-five  combinations  of  errors  were  calculated. 
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The  standard  deviations  and  maximum  errors  were  incorporated  into 
plots  of  the  turbulence  parameters.  A plot  was  made  for  each  parameter 
for  each  latitude  for  each  of  three  seasons  (winter,  summer,  annual),  with 
the  corresponding  standard  deviation  at  ten  kilometer  intervals.  For 
the  plots  of  the  Richardson  number,  the  standard  deviation  was  not 
necessarily  used.  Since  the  maximum  possible  error  was  to  be  displayed, 
the  greater  of  the  standard  deviations  or  the  experimental  errors  was 
used.  Plots  were  also  made  of  the  worst  possible  case  for  each  latitude. 

After  the  calculation  of  the  turbulence  parameters  had  been  done 
on  a summer,  winter  and  annual  basis,  further  breakdown  was  done  for  one 
location.  There  was  enough  data  from  Wallops  Island  to  have  statistically 
significant  results  if  done  seasonally.  The  data  was  run  for  four  time 
periods  corresponding  to  four  seasons. 

Additional  data  was  obtained  in  an  effort  to  verify  the  results 
obtained  for  the  four  latitudes.  Another  set  of  data  from  a latitude 
close  to  that  of  Wallops  Island  was  used.  However,  the  results  were  not 
as  expected.  Originally,  the  discrepancy  was  attributed  to  the  presence 
of  high  mountains  near  the  new  site  which  could  have  created  disturbances. 
Consultation  with  the  originator  of  the  data  revealed  that  two  differ- 
ent models  were  used  to  reduce  the  data,  one  for  higher  altitudes 
and  one  for  lower  altitudes.  This  data  was  discarded  and  additional  data 
requested.  This  sequence  is  reported  to  indicate  that  there  have  been 
problems  in  the  past  with  the  data  base.  Care  must  be  taken  in  the  inter- 
pretation of  results  obtained  from  the  data  available. 

Another  comparison  was  attempted  with  additional  data,  predom- 
inantly low  altitude.  This  data  was  available  for  altitudes  between  twenty 
and  sixty  kilometers  for  two  locations.  One  location  was  Wallops  Island. 
Hence  there  was  an  overlap  with  existing  altitude  data  from  that  site. 

The  other  location  was  Ascension  Island  which  has  a latitude  similar  to 
that  at  Natal.  (Data  from  Natal  had  been  previously  processed  and 
analyzed.)  When  these  two  data  sets  were  processed,  no  occurrences  of 
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turbulence  were  found  despite  the  fact  that  some  instances  appeared  at 
the  low  end  of  the  high  altitude  data.  There  are  two  possible  explana- 
tions for  these  results.  First,  and  most  likely,  is  that  since  this  data 
was  smoothed,  the  smoothing  process  masked  or  eliminated  the  data  which 
would  have  resulted  in  turbulence.  The  second  possibility  is  that  the 
higher  altitude  data  was  incorrect  since  it  was  at  the  low  end  of  its 
range.  Subsequent  investigation  indicates  that  the  smoothing  done  was 
minimal  and  that  is  probably  did  not  significantly  affect  the  data. 

Additional  plots  were  requested  and  generated  for  this  task.  The 
previously  plotted  turbulence  parameters  were  put  through  a modified  ver- 
sion of  an  existing  interpolation  routine  to  provide  a smoothed  least 
squares  fit  of  the  data.  These  plots  also  displayed  the  weighting  factor, 
or  number  of  data  sets  on  which  the  point  was  based,  next  to  the  point. 
Further,  plots  of  the  wind  components  were  made.  There  were  several  methods 
used  in  displaying  the  wind  components. 

The  first  set  of  plots  showed  the  wind  components  for  all  experi- 
mental shots  plotted  together.  These  showed  what  was  anticipated:  namely 

that  tne  east-west  component  exhibited  a fairly  constant  value  over  all 
altitudes;  whereas,  the  range  of  the  north-south  component  increased  with 
increased  altitude.  All  latitudes  showed  this  phenomenon,  with  the  results 
from  Natal  being  the  most  pronounced. 

Other  plots  were  made  of  the  wind  components.  Each  shot  was 
plotted  alone  with  a ten  point  running  average.  Then,  the  average  was 
subtracted  from  the  component  and  the  result  plotted.  A program  was 
also  written  to  produce  histograms  of  the  wind  components.  This  was 
done  to  determine  whether  the  distribution  of  the  wind  follows  a normal 
distribution.  Further  processing  and  analysis  is  anticipated  on  this  task. 

Mary  of  the  results  and  plots  from  this  effort  were  presented  at 
o NA'iV  .onference  r April  1977.  Additional  results  and  plots  will  be 
published  shortly. 


Additional  high  and  low  altitude  data  for  1969-1975  has  been 
received  and  is  being  processed.  The  high  altitude  tapes  have  been 
read  and  processing  will  be  performed.  The  results  of  this  initial 
processing  will  determine  what  further  processing  will  be  done.  The 
preliminary  reading  of  the  low  altitude  data  is  currently  underway. 


U6 


25. 


Ionospheric  Mapping 


INITIATOR:  C.  RUSH 

PROBLEM  NO:  4927  PROJECT  NO:  7663 


BACKGROUND 

The  objective  of  this  effort  is  to  determine  the  accuracy  of 
current  recommended  techniques  for  forecasting  and  specifying  ionospheric 
parameters.  Predicted  fQF2  values  will  be  compared  with  observed  data. 
Measured  monthly  median  (hourly)  ionospheric  parameters  are  available  from 
various  locations  around  the  world.  Tapes  of  these  data  have  been  pre- 
pared for  analysis.  Once  the  accuracy  of  the  world-standard  has  been 
evaluated  and  described  as  a function  of  local  time,  season,  solar  cycle, 
and  geographical  location,  it  will  be  possible  to  improve  upon  the  fore- 
casting accuracy  of  the  original  computations  by  incorporating  corrections 
based  on  the  measured  data. 

This  task  is  not  completed.  This  following  is  a summary  of  what 
has  been  accomplished  and  what  additional  work  is  needed  for  the  future 
analysis . 


ANALYSIS  AND  RESULTS 


Specific  tasks  being  undertaken  are: 

Task  1.  Sort  each  tape  according  to  month  (original  tape 
arranged  according  to  station) . Retain  only  the 
values  associated  with  the  fQF2  values.  Four  tapes 
have  been  processed,  and  two  additional  tapes  have 
been  submitted  for  processing.  During  preparation 
of  the  data  tapes  for  processing,  numerous  errors 
have  been  detected  on  the  tape.  Data  records  are 
missing.  To  provide  for  a "clean  data  tape"  for 
the  analysis  program,  dummy  records  were  added. 
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Header  records  are  missing.  For  this  case  the 
site,  its  coordinates,  etc.,  are  being  obtained 
from  tables  supplied;  and  this  information  is 
being  added  to  the  data  on  the  tape.  This  infor- 
mation, header  and  data  records,  have  been  added 
to  the  four  ta'pes  already  sorted.  The  other  two 
tapes  will  be  supplemented  with  the  information 
and  processed. 

Further,  problems  have  occurred  with  the  reading 
of  one  of  the  data  tapes  supplied.  The  source 
ot  the  data,  National  Bureau  of  Standards,  was 
;onta  ted.  I t became  apparent  that  the  tape 
■ nt  was  mislabeled.  NBS  agreed  to  send  another 
tape  replacing  the  original  one. 

Task  2.  Provide  a listing,  by  month,  of  stations  from 
which  data  are  available.  The  four  tapes  men- 
tioned previously  have  been  listed.  Data  has 
been  listed  for  the  years,  1970,  1972,  1973,  and 
1974.  After  the  preparation  of  the  data  tapes 
additional  data  will  be  available  for  1966  and 
1971 . 


Task  3.  For  each  month  (and  appropriate  sunspot  number) 

compute  the  difference,  for  each  hour,  between  the 
observed  median  f0l'2  ant*  that  predicted  using  the 
■oeffi  ■lent:  contained  in  the  modified  "ITS  Model" 

that-  is  :urrently  being  used  in  the  ISIS  2 analysis 
Record  the  diffor"nce  between  observation  and  pre- 
diction, together  with  the  number  of  observations 
that  went  into  the  computation  of  the  median. 
Software  has  been  developed  to  perform  this  task. 
Again  four  tapes  have  been  prepared  for  future 
analysis. 

Task  4.  Arrange  the  output  on  a tape  so  that  the  following 
can  easily  be  obtained: 

a.  Plots  of  differences  for  a given  station- 
month  as  a function  of  time  - 24  values 
for  each  of  the  24  hours. 

b.  Values  of  differences  (or  differences  squared) 
placed  at  the  appropriate  station  coordinates 
(or  geomagnetic  coordinates)  for  each  hour  of 
local  mean  or  universal  time.  This  is  equiva- 
lent to  placing  individual  values  of  the 
differences  at  the  appropriate  locations  on 

a global  map. 

e.  RMS  values  of  differences  computed  from  spe- 
cified grouping(s)  of  stations. 
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Figure  25-1  illustrates  the  geographical  distribution  of  the 
stations  supplying  data  for  a particular  month.  This  plot  shows  that 
data  is  available  from  locations  around  the  world.  (It  was  not  intended 
to  be  accurate  within  each  15°  x 15°  block.  An  "X"  does  not  indicate 
exact  location.  Rather  an  X indicates  a site  is  located  within  the  block, 
two  X's  indicate  two  sites  within  the  block.)  After  analyzing  the  dis- 
tribution of  sites  it  was  decided: 

1.  Not  to  perform  a computer  contour  interpola- 
tion, at  present.  This  is  because  there  is 
no  information  for  ocean  areas. 

2.  To  determine  relationships  between  sites 
within  blocks.  That  is,  determine  if 
differences  between  predicted  and  observed 
f0F2  for  sites  relatively  close  to  each 
other  can  be  modeled  by  using  only  one  site. 


The  analysis  will  continue  by  developing  software  to  plot  on  a Mercator 
projection  at  the  site  locations: 


1. 


The  difference  between  predicted 
f0F2  for  all  sites  for  the  hours 
local  time  and 


and  observed 
00,  03,  06,  ... 


The  difference  between  predicted  and  observed 
f0F2  for  all  sites  for  the  hours  00,  03,  06,  ... 
universal  time. 


A plot  will  be  produced  for  each  hour/month  combination.  Characters  will 
be  displayed  to  indicate  a range  of  differences.  This  effort  is  a sub- 
task of  Task  4 described  previously.  Table  25-1  is  representative  of  the 
data,  both  predicted  and  observed. 
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26.  Correlation  Studies  For  The  Prediction 
Of  Scintillation  And  Total  Flectron  Content 

INITIATOR:  J.  AARONS 

PROBLEM  NO.:  4918  PROJECT  NO.:  4643 


BACKGROUND 


ASEC  has  been  involved  with  a ta-k  whose  purpose  is  to  determine 
if  scintillation  and  total  electron  content  'TEC)  at  a station  can  be  pre- 
dicted from  the  scintillation  and/or  TEC  at  another  station.  The  aim  is  to 
find  forecasting  tools  for  AWS  in  problems  of  radar  corrections  and  in 
irregularity  formation.  Specifically,  the  first  subtask  is  to  develop  a 
series  of  programs  to  calculate  correlation  functions  in  support  of  scin- 
tillation studies  for  AFSATCOM  and  total  electron  content  studies  in  support 
of  ADCOM  and  GPS. 

The  material  for  the  studies  consists  of  series  of  values  of 
scintillation  indices  and  of  total  electron  content  values.  The  volume 
of  data  extends  over  several  years  and  consists  of  measurements  for  many 
satellites  and  for  several  locations.  Basically,  correlation  functions 
(auto  and  cross)  for  these  15-minute  indices  are  required  as  a first  approach 
in  the  analysis. 

Previous  studies  have  addressed  and  documented  a high  latitude 
model  of  scintillation  excursion  based  on  observations  of  the  scintilla- 
tions of  beacons  from  synchronous  satellites.^  Equations  have  been 
developed  to  predict  scintillations  as  a function  of  local  time,  magnetic 
index,  solar  flux,  and  month  of  the  year.  Other  studies  have  shown  clear 
seasonal  maximum  of  the  occurrence  of  equatorial  scintillations.  Diurnal 


"A  High-Latitude  Empirical  Model  of  Scintillation  Excursions: 
Phase  I",  Jules  Aarons  et  al,  17  September  1976,  AFGL-TR-0310. 
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occurrence  of  scintillation  peaks  before  midnight  for  both  quiet  and 

2 

magnetically  disturbed  days.  Based  upon  these  studies  and  others,  it 
appeared  possible  to  predict  scintillation  and  total  electron  content  at 
a station,  given  knowledge  of  scintillation  and/or  TEC  at  another,  spatially 
separated  station.  Consequently  this  prediction  task  was  undertaken. 


ANALYSIS  AND  RESULTS 

The  initial  objective  was  to  develop  a software  package  to  perform 
correlation  calculations.  The  material  for  the  studies  consists  of  series 
of  values  of  scintillation  indices  and  of  total  electron  content  values. 

The  volume  of  data  extends  over  several  years  and  consists  of  measurements 
for  many  satellites  and  for  several  locations.  The  data  consist  of  scin- 
tillation values  in  dB  or  in  scintillation  index.  They  consist  of  TEC 
absolute  values.  Basically,  correlation  functions  (auto  and  cross)  for 
these  15-minute  indices  are  required. 


The  following  is  a summary  of  work  performed  on  this  problem  to 
date.  Although  the  indices  are  of  15  minutes  of  data,  further  calculations 
may  involve  one  hour  mean  indices  for  auto  and  cross  correlations.  Provi- 
sions for  this  data  handling  have  been  included  in  the  software.  It  is 
presently  possible  to  obtain  nightly  values  and  use  these  values  to  com- 
pute auto  and  crosscorrelation  functions.  If  desired,  the  four  highest 
values  for  the  period  2200-0400  UT  may  be  specified  and  used  as  the  basis 
for  the  auto  or  crosscorrelation  functions.  It  is  possible  to  soecify 
auto  and  crosscorrelation  functions  for  one  season  or  one  month.  Careful 
distinction  is  made  in  the  data  between  no  data  and  zero  values.  (Portions 
of  the  data  supplied  were  supplemented  with  a code  to  identify  missing  data.) 


2 

"Equatorial  Scintillations: 
1976,  AFGL-TR-76-0078 . 


A Review",  Jules  Aarons,  13  April 
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Software  has  been  jeie- rated  to  alculate  auto  and  crossrorrela- 
tion  functions.  Both  five-day  lag  and  lb-hour  lag  plots  are  being  pro- 
duced for  future  analysis.  s.imi  1<-  plots  are  shown  in  figures  26-1,  26-2, 
26-3,  and  26-4.  Total  data,  and  monthly  and  seasonal  auto  and  cross- 
correlation functions  are  being  produced  for  three  sites.  They  are  Narssar- 
ssuaq.  Goose  Bay,  and  Sagamore  Hill.  In  the  calculations,  the  maqnetic 
index  can  be  quantized  in  steps  of  selectable  size.  Initially,  two  scales 
of  granularity  were  used: 

a)  A division  with  Kp  either  in  the  range  0 to  3+  or 
in  the  range  of  4-  to  9 

b)  No  division  with  respect  to  Kp. 


Using  only  nightly  data  has  produced  results  which  show  diurnal  variations 
about  midnight,  which  was  to  be  expected. 


One  of  the  difficulties  has  been  the  handling  of  large  amounts  of 
data  for  multi-seasonal  and  total  time  auto  and  crosscorrelations.  The 
problem  has  been  resolved  with  the  inclusion  of  an  algorithm  which  combines 
means,  varriances,  auto  covariances  or  cross-covariances  from  the  individual 
time  groups  of  data.  For  disjoint  data  groups,  the  equations  used  for  the 
combinations  of  means,  variances,  auto-covariance  or  total  cross-covariance 
are : 


Total  Mean:  X 

X. 

l 

N. 

1 

N 

Total  Variance: 


EN.X. 
l l 

N 

= mean  of  individual  group 

= number  of  occurrences  in  individual  group 

= LNi  2 - - 2 

, E (N.o.  ) + E N.  (X.  - X) 

^ 11  11 
a = — 


N 


where  o . 


variance  of  individual  group. 
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Total  auto-covariance:  TACV . = |E  N.ACV. ,+  £[M. (X.  - X)^"l  + (X,  - X) 
3 L 1 ij  L l l J l 

[ (E( 


where  = number  of  occurrences  in  group  i 

X^  = means  of  group  i 
X = overall  mean 

NT  = number  of  terms  summed  in  the  group  i 
ACVij  = auto  covariance  for  groups  i,  lag  j 


Hence  the  total  auto  correlation  is  determined  by 
TACV. 

AC.  = 

3 


A similar  computation  for  the  crosscorrelation  becomes: 


Total  cross  covariance: 

TCCV . = |e  N.CCV.  . + £ [m.  (X.  - X^)  (Y.  - V 1 

3 L i3  Mi  1 i _ Mi 

+ E (X . - X)  E (Y.  - Y.)  + E ( Y . - Y)  E (X.  - X)  /£n. 

l ,,li  l ,l  I l 

1=1  1=1  J 

for  groups  Y^,  Y^ , Y^  lagging  groups  X^,  X^ 

N.  = number  of  occurrences  used  to  calculate  cross- 
covariance for  group  Y^  lagging  X^ 

X^  = means  of  group  X^ 

Y.  = means  of  group  Y^ 

CCVij  = cross-covariances  for  group  i (Y^  lagging  XJ , lag  j 
M.  = number  of  terms  used  in  the  calculation 

l 

X = overall  mean  of  the  X^  's 
Y = overall  mean  of  the  Y^  's 


X..  ) ) - 2X 


iMJj' 
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Hence, the  crosscorrelation  is  determined  by 


CC 


TCCV  _ 


OX 


OY 


13 

2 


r 2 2 -1 1/2 

[aX  °Y  ] 
total  variance  for  group  X 

total  variance  for  group  Y 


In  addition  to  these  algorithms,  the  individual  group  (season) 
auto-covar lances , number  of  data  points,  means,  etc.,  are  saved  for  future 
use.  Hence,  this  method  facilitates  the  inclusion  of  additional  data  when 
it  becomes  available.  Past  results  are  not  recalculated;  rather,  they  are 
new  terms  or  constants  which  are  added  to  the  equations.  Overall  data 
handling  is  facilitated,  and  computer  time  is  saved.  Modifications  to  the 
software  to  include  data  from  overlapping  groups  have  been  formulated  and 
■uld  be  included. 


The  first  phase  of  the  analysis  software  has  been  completed. 

This  software  has  been  used  to  produce  auto  and  crosscorrelations  based 
on  the  scintillation  index  at  the  three  sites.  At  present,  some  3 00  sets 

of  autocorrelation  plots  and  100  sets  of  crosscorrelation  plots  have 
been  produced  for  future  analysis. 

Future  work  will  include  the  completion  of  the  calculation  of  all 
auto  and  crosscorrelation  functions  for  the  TEC  and  crosscorrelat ions 

between  TEC  and  scintillation  index  for  all  combinations  (i.e.,  seasonal, 
nightly,  yearly).  Following  this,  will  be  the  calculation  of  monthly 
means  and  the  generation  of  correlation  functions  of  the  deviations  from 
monthly  means  and  the  crosscorrelations  for  deviations  at  several  sta- 
tion:.. Most  of  the  software  necessary  to  calculate  the  correlation  func- 
tions of  the  deviations  from  monthly  means  is  available.  Additions  to 
'■xi  .ting  software  in  use  will  be  necessary  to  complete  the  software  package. 
Upon  completion  of  the  correlation  function  calculations  and  analysis 
thereof,  equations  will  be  formulated  to  predict  the  scintillation  index 
and  total  electron  content  at  a specified  geographic  location  given  values 
of  scint illation  and/or  TEC  at  another  station. 
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Figure  26-2  Goose  Bay  Autocorrelation  (Hours) 


27.  Aerosol  Spectroscopy 


INITIATOR:  F.  VOLZ 

PROBLEM  NO. : 4928  PROJECT  NO. : 7621 


BACKGROUND 

Particle  counts  are  being  measured  by  the  AFGL/OPA  group  at  Hanscom 
AFB.  It  is  desired  that  the  particle  counts  be  obtained  with  respect  to 
size  over  a period  of  time.  This  data  is  obtained  utilizing  a DAS-64  Dual 
Differential  Analog  Module,  and  values  are  recorded  on  a 9-track  tape. 

The  particles  are  being  counted  in  six  different  ranges.  These 
>anges,  in  turn,  are  subdivided  into  15  equal  groups,  referred  to  as  bins. 
The  particle  counts  are  recorded  at  a sampling  period  determined  by  the 
operator.  This  period  is  between  .01  seconds  to  1 second.  Along  with 
the  cour.' s for  the  various  ranges,  the  date,  time,  counting  rate,  and 
range  size  are  stored  on  the  tape  for  reference  and  future  analysis. 

Data  is  obtained  either  in  automatic  mode,  where  all  ranges  are 
sample  sequentially,  or  manually,  where  a particular  ranae  may  be  sampled. 
All  analysis  here  assumes  automatic  mode. 


ANALYSIS  AND  RESULTS 

The  data  as  it  appears  on  the  tape  is  described  in  Figure  27-1. 
Four  bits  of  each  word  represent  a field.  Hence  it  was  necessary  to  write 
a program  to  unpack  the  data  on  the  tape  into  CDC  6600  words  containing 
4 bits  of  information.  These  fields  were  then  combined  using  the  scale 
factor  described  in  Figure  27-1  to  obtain  the  data  values. 
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Figure  27-1  DAS-64  Tape  Format 
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Software  has  been  written  to  read  a 7 inch  phase  encoded  9-track 
hexadecimal  tape  with  aerosol  count  data  and  to  rearrange  and  display  this 
data.  The  data  for  log-log  plots  of  counts  versus  size  of  bins  has  been 
normalized.  In  addition  to  providing  the  counts  for  the  fifteen  different 
size  bins  for  each  of  six  measuring  ranges,  the  software  performs  the  fol- 
lowing operations. 

1.  Calculation  of  absolute  particle  concentration 

2.  Determination  of  bin  size  distribution 

3.  Definition  of  bin  size 

4.  Smoothing  of  data  to  extract  large  statistical 
fluctuation 

5.  Printing  of  information  (format  described  in 
Figure  27-2  and  27-3) 

6.  Plotting  of  information  (Figure  27-4) 

7.  Provision  for  count  data  to  be  input  via  cards. 


The  following  is  a discussion  of  the  equations  used  and  the 
output  obtained  from  the  software.  Particle  counts  in  a certain 
size  bin  (i  = 1 to  15)  of  one  of  the  6 measuring  ranges  are  related 
to  absolute  particle  concentration  dir  by 


dN 

l 


C. 

1 

T • SVR 


(cm 


) 


C = number  of  counts 


T = Collection  duration  (seconds) 


SVR  = Sampling  volume  rate  (cm  sec  ) 


This  concentration  refers  to  the  size  (diameter)  interval 


dD.  = D.  , - D.  or  radius  interval  dR . = dD./2 
l l + 1 l li 


In  a flat  size  distribution,  or  narrow  bins,  the. median  radius  R . is 

mi 

calculated  with  the  equation  R . = (D.  , + D.)/4.  However,  where  the 

mi  l + 1 l 
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bin  size  is  wide  and  the  size  distribution  steep,  an  effective  radius 
may  be  defined  as: 


R 

ei 


[D.  + K (D . ' 

L l l + 


1 


D.)]/2 


where  input  parameter  K (0  < K < 1);  in  natural  haze,  K ~ 0.2 


R = R for  K = 0.5 
m e 


The  range  sizes  are  defined  below: 

Range  Size  in  microns  (y)  Bin  Size  in  U (Bin  #1  thru  15) 


A3 

0.15  to 

0.3 

(.01) 

A2 

0.2 

to 

0.65 

(.03) 

Al 

0.35  to 

1.1 

(.05) 

AO 

0.6 

to 

3.0 

(.16) 

Bl 

0.6 

to 

3.0 

(.16) 

BO 

1.0 

to 

20.0 

(1.27) 

The  effective  radius  is  calculated  for  Ranges  AO,  Bl,  BO,  for  bin 
sizes  1 through  6.  All  other  ranges  use  the  median  radius. 


It  was  desired  to  calculate  the  size  distributions  in  terms  of 

dN  /dlogr  . 
l l 


dlogr.  = 0.434  (£n  R . . -£n  R .) 

^ i ei  + l ei 


2.3  C. 

Therefore  dN./dlogr.  = --  - -f  _ -;n  R;}  't'T  SVR 


Range  cind  size  data  will  be  different  for  various  instruments.  This  infor- 
mation is  input  via  data  cards. 
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Low  counts  show  large 

statistical  fluctuations. 

During  conversion  to 

size 

distributions , 

they  are  smoothed  as 

follows : 

SC.  = 

X 

(C. 

l - 

, + 4C.  + C. 

1 l l + 

l)/6 

Present  count 

limits 

for  smoothing 

have  been  set 

for  values  < 20. 

Counts  for  particular  range  and  bin  sizes  are  summed  for 

periods  of  time  determined  on  input 

cards,  and  this  is  output  to  the 

printer  as  in 

Figure 

27-2. 

DATE  TIME 

RANGE 

BIN 

0 1 

2 3 

. . 15 

3-4 

Counts 

8500  3000 

980 

2-4 

2210 

Figure  27-2 

Count  Data 

An  option  provides  count  data  from 

input  cards  if 

desired. 

Count  listing 

is  followed  by  a dN/dlogr  listing. 

An  example  is  given 

in 

Figure  27-3: 

DATE  TIME 

RANGE 

BIN 

0 

1 2 

3 

171830 

1-4 

Radius 

.015 

.02  .025 

dN/dlogr 

3.32  + 2 

3.0  + 1 

171900 

0-4 

Radius 

.65 

.77 

dN/dlogr 

2.2  + 1 

Figure  27-3 

Listing  of  dN/dlogr 

Log  - log  plots  of  the  dN/dlogR  versus  radius  are  obtained  with 
all  ranges  on  one  plot.  Each  range  identified  by  a separate  symbol.  The 
method  of  processing  of  data  taken  in  a manual  mode  will  be  decided  at  a future 
date . 
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Calculations  and  processing  have  been  made  for  sample  data. 

Further  modification  to  the  processing  software  may  be  necessary.  Smoothing 
of  output  to  be  displayed  has  been  discussed.  Algorithms  necessary  to  inter- 
polate output  are  available  and  may  be  incorporated.  Certain  anomalies  have 
been  detected  in  the  instruments.  Once  these  have  been  corrected  by  the 
vendor,  interpolation  may  be  applied  to  the  output  before  plotting.  This 
program  will  be  used  for  the  comparison  of  data  obtained  from  NATO  aerosol 
spectrometers.  Figure  27-4  is  a sample  of  the  output  plot  being  generated. 


iqure  27-4  Cutout  from  Aerosol  Spectrometer 


28.  German  Surface  Weather  Data 


INITIATOR:  I.  GRINGORTEN 


PROBLEM  NO:  4931  PROJECT  NO:  8624 


BACKGROUND 

Twenty  seven  magnetic  tapes  containing  weather  information  from 
26  German  weather  stations  have  been  submitted  for  processing  and  analysis. 

It  is  desired  to  extract  data  such  as  ceiling,  visibility,  date,  time, 
station  coordinates  and  station  numbers  from  these  tapes,  and  to  store 
this  information,  reformatted,  on  separate  tapes.  This  data  will  be  used 
to  calculate  probabilities  distributions  for  visibility  and  ceiling.  These 
tapes  were  generated  on  a 32  bit  machine  and  each  field  is  composed  of 
either  16  or  32  bits  of  information.  A conversion  routine  unpacks 
60  bits  of  input  information  into  a 60  bit  word  containing  16  bits  of  in- 
formation. In  the  way  each  field  could  be  easily  represented  by  either 
one  or  two  words  of  information.  Since  all  data  was  in  integer  form,  unless 
otherwise  stated,  there  was  no  need  for  floating  point  conversion.  Since 
four  CDC  6600  sixty-bit  words  equal  fifteen  16-bit  words,  a loop  was 
set  up  in  the  program  to  convert  four  CDC  input  words  of  information  into 
15  words  containing  16  bits  of  information  with  zero-fill  on  the  left. 

ANALYSIS  AND  RESULTS 

Problems  arose  with  the  control  section  of  the  data.  According 
to  the  data  documentation, the  variable  in  field  CO 1,  record  length,  was  to 
contain  the  total  number  of  bytes  in  the  record.  However,  this  was  not 
found  to  be  the  case.  Since  the  block  size  and  record  size  is  a variable, 
it  was  necessary  to  search  for  the  records  by  keying  on  an  easily  iden- 
tifiable parameter.  Since  the  first  record  of  each  block  always  had  the 
same  length,  except  for  the  "additional  Data  Section,'  which  appears  after 
the  needed  data,  the  year,  month,  day,  hour  could  be  obtained  from  the 
first  record  of  each  block.  By  searching  for  these  parameters  in  the 
subsequent  record,  all  records  in  each  block  could  be  obtained. 


The  latitude  and  longitude  were  converted  from  degrees  and 
minutes  to  degrees  to  nearest  hundredth  by  dividing  the  minutes  by  60 
and  adding  this  fraction  to  the  number  of  degrees.  The  visibility  value  on 
tape  in  meters  was  divided  by  402,  and  if  it  exceeds  100,  it  was  set  to  100. 
The  value  402  meters  is  the  largest  expected  value. 


Multi-file  tapes  were  set  up  for  the  data  extracted  from  the 
tapes, converted  and  stored  in  the  following  order. 


ORDER 

FIELD 

DESCRIPTION 

CONVERT  TO 

SYMBOL 

1 

C06 

Month 

Same  (2  digits) 

MONTH 

2 

C09 

Hour 

01  to  24 
(2  digits) 

NHOUR 

3 

C04 

Block-station 

Numbered  1 to  26 

1ST 

4 

C12 

Latitude 

Degrees  to 
nearest  100th 

PLAT (1ST) 

5 

C13 

Longitude 

Degrees  to 
nearest  100th 

PLONG (1ST) 

6 

CO  5 

Year 

Last  two  digits 
of  year 

IYEAR 

7 

C07 

Day 

01  to  31 

MDAY 

8 

M30 

Ceiling 

Same  WMO  code 

IH 

9 

M12 

Visibility  (V) 

v = Integer  of 

(V/402)  for  V < 40200 
= 100  for  V ■>  40200 

The  above  data  will  be  used  along(  additional  software  to  perform 
the  following: 


a.  Calculate  the  probability  distributions  of 
cloud  ceiling  heights  (>  H)  and  visibility 
(>_V)  at  each  of  the  26  stations  and  the 
equivalent  normal  deviates  (e.n.d.)  for  the 
given  heights  and  visibilities. 
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b.  Calculate  the  correlation  coefficients  (c.c.) 
for  all  combinations  of  the  e.n.d.'s  of 
ceiling  (cig)  and  visibility  (vis)  at  the 

26  stations. 

c.  Calculate  the  joint  probability  of  cig  and 
vis  at  each  station. 

d.  Determine  the  values  of  the  parameter,  known 
as  scale  distance,  for  cig  and  for  vis,  in 
the  model  formula  for  estimating  c.c. 

e.  Calculate  the  conditional  probability  of 
cig  and  vis  and  the  ]oint  probability  at 
a point  (4>,  A)  . 

f.  Determine  the  probability  distribution  of 
the  minimum  of  cig  and  vis  in  areas  of 
varying  size. 


Portions  of  the  original  27  tapes  were  copied  onto  other  in-house 
tapes  for  future  reference.  Additional  information  may  be  extracted  at  a 
future  date. 
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29.  Floating  Point  Arithmetic  For  CDC  6600 


Floating  point  arithmetic  in  the  CDC  6600  takes  advantage  of  the 
ability  to  express  a number  with  the  general  expression  kB*1,  where 

k = coefficient 
B = base  number, 


and 


n = exponent,  or  power  to  which  the  base  number  is  raised. 


The  base  number  is  constant  for  binary-coded  quantities  and 
is  not  included  in  the  general  format.  The  60-bit  floating-point  format 
is  shown  in  Figure  29-1  .The  binary  point  is  considered  to  be  to  the  right 
of  the  coefficient,  thereby  providing  a 48-bit  integer  coefficient,  the 
equivalent  of  approximately  14  decimal  digits.  The  sign  of  the  coefficient 


is  carried  in  the  highest-order  bit  of  the  packed  word.  Negative  numbers 
are  represented  in  one's  complement  notation. 

8 0 0 7 


COEFFICIENT  BIASED  INTEGER 

SIGN  EXPONENT  COEFFICIENT 


59  58  48  47 


{ I 

0 

BINARY 

POINT 


ricmre  29-1  60-bit  Floating  Point  Format 

The  11-bit  exponent  carries  a bias  of  210  (2000g)  when  packed  in 
the  floating  point  word  (biased  exponent  sometimes  referred  to  as  character- 
istic) . The  bias  is  removed  when  the  word  is  unpacked  for  computation  and 
restored  when  a word  is  packed  into  floating  format.  For  complete  range 
of  permissible  values  see  the  CDC  6600  SCOPE  reference  manual. 
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